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Abstract. We perform direct numerical simulations of the Hamiltonian mean field (HMF) model starting 
from non-magnetized initial conditions with a velocity distribution that is (i) gaussian, (ii) semi-elliptical, 
and (iii) waterbag. Below a critical energy Ec, depending on the initial condition, this distribution is 
Vlasov dynamically unstable. The system undergoes a process of violent relaxation and quickly reaches a 
quasi-stationary state (QSS). We find that the distribution function of this QSS can be conveniently fitted 
by a polytrope with index (i) n = 2, (ii) n = 1, and (iii) n = 1/2. Using the values of these indices, we are 
able to determine the physical caloric curve Tkin{E) and explain the negative kinetic specific heat region 
Ckin = dE/dTkin < observed in the numerical simulations. At low energies, we find that the system 
takes a "core-halo" structure. The core corresponds to the pure polytrope discussed above but it is now 
surrounded by a halo of particles. In case (iii), we recover the "uniform" core-halo structure previously found 
by Pakter & Levin [Phys. Rev. Lett. 106, 200603 (2011)]. We also consider unsteady initial conditions with 
magnetization Mo = 1 and isotropic waterbag distribution and report the complex dynamics of the system 
creating phase space holes and dense filaments. We show that the kinetic caloric curve is approximately 
constant, corresponding to a polytrope with index no — 3.56 (we also mention the presence of an unexpected 
hump). Finally, we consider the collisional evolution of an initially Vlasov stable distribution, and show 
that the time-evolving distribution function f{v,t) can be fitted by a sequence of polytropic distributions 
with a time-dependent index n{t) both in the non-magnetized and magnetized regimes. These numerical 
results show that polytropic distributions (also called Tsallis distributions) provide in many cases a good 
fit of the QSSs. They may even be the rule rather than the exception. However, in order to moderate our 
message, we also report a case where the Lynden-Bell theory (which assumes ergodicity or efficient mixing) 
provides an excellent prediction of an inhomogeneous QSS. We therefore conclude that both Lynden-Bell 
and Tsallis distributions may be useful to describe QSSs depending on the efficiency of mixing. 

PACS. 5.20.-y Classical statistical mechanics - 05. 45. -a Nonlinear dynamics and chaos - 05.20.Dd Kinetic 
theory - 64.60.De Statistical mechanics of model systems 



1 Introduction 

Systems with long-range interactions are numerous in 
nature. They include, for example, self-gravitating sys- 
tems and two-dimensional (2D) vortices which are sys- 
tems of considerable interest |l]. These systems may 
be trapped in long-lasting non-equilibrium states, called 
quasi-stationary states (QSS), whose lifetime diverges 
with the number of particles N. These QSSs correspond 
to galaxies in astrophysics and large-scale vortices (e.g. 
Jupiter's great red spot) in 2D hydrodynamics. There- 
fore, in many cases of physical interest, the system does 
not reach the Boltzmann distribution but remains stuck 
in a non-Boltzmannian QSS. This is the case of elliptical 
galaxies in stellar dynamics because the collisional relax- 
ation time exceeds the age of the universe by many orders 
of magnitude. This is also the case in 2D geophysical and 



astrophysical flows because the viscous time is generally 
much larger than the turnover time of a large-scale vortex. 
These QSSs are known to be stable steady states of the 
Vlasov (or 2D Euler) equation on a coarse-grained scale. 
The Vlasov equation describes the "coUisionless" evolu- 
tion of the system before "collisions" (more precisely cor- 
relations, finite N effects, granularities...) drive the sys- 
tem towards Boltzmann's statistical equilibrium. Since 
the Vlasov equation admits an infinite number of steady 
states, the prediction of the QSS that is actually selected 
by the system is difhcult. Our understanding of these QSSs 
is still incomplete. 

A toy model of systems with long-range interactions, 
called the Hamiltonian mean field (HMF) model, has been 
actively studied in statistical mechanics [2J. It consists of 
N particles moving on a ring and interacting via a cosine 
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potentiaQ At statistical equilibrium, this system displays 
a second order phase transition between a spatially homo- 
geneous (non-magnetized A/ = 0) phase and a spatially 
inhomogeneous (magnetized M 7^ 0) phase. The magne- 
tized phase appears below the critical energy Ec = 3/4 
or below the critical temperature Tc — 1/2. Antoni & 
RufFo 8 carried out direct numerical simulations of the 
HMF model. They started from an initial condition in 
which all the particles are located at 6* = (correspond- 
ing to a magnetization Mq = 1) with a waterbag velocity 
distribution. They determined the physical caloric curve 
Tkin{E) giving the average kinetic temperature as a func- 
tion of the energy. They compared their numerical results 
to the theoretical caloric curve T(E) corresponding to the 
Boltzmann equilibrium and reported several "anomalies" . 
In particular, the phase transition takes place at an en- 
ergy sensibly smaller than Ec = 3/4 and the numerical 
caloric curve Tkin{E) presents a region of negative kinetic 
specific heats, unlike the theoretical caloric curve T{E) 
corresponding to the Boltzmann equilibrium. They under- 
stood that these discrepancies are due to the fact that the 
observed structures are out-of-equilibrium QSSs. These re- 
sults were confirmed by Latora et al. [ll|12j who showed 
that the lifetime of these QSSs diverges with and that 
their distribution functions are non-Boltzmannian. They 
also observed many other anomalies in the region of neg- 
ative kinetic specific heats such as anomalous diffusion. 
Levy walks, aging, and dynamical correlations in phase- 
space. The observation of non-Boltzmannian QSSs was a 
surprise in the community of statistical mechanics. La- 
tora et al. |11I12| proposed to interpret these QSSs in 
terms of Tsallis generalized thermodynamics ^13 . In par- 
ticular, they tried to fit the QSS a,t E = 0.69 by a q- 
distribution with a power-law tail. To make the distribu- 
tion normalizable, they introduced a cut-off at large veloc- 
ities. While their study definitely shows that the QSS is 
non-Boltzmannian, their procedure is not very convincing 
and their fit is relatively poor. 

The situation changed after the conference in Les 
Houches in 2002 where it was indicated ^4] that QSSs 
were previously observed in stellar dynamics and 2D tur- 
bulence. In these domains, the QSSs are interpreted in 
terms of I/ynden-Bell's statistical theory of violent relax- 
ation [15 Fl The Lyndcn-Bcll theory determines the statis- 
tical equilibrium state of the Vlasov equation, taking into 
account all the constraints of the "coUisionless" dynam- 
ics, in particular the conservation of the Casimirs. If the 
system is ergodic, the QSS coincides with the Lynden-Bell 
distribution (most probable state). In the two-levels case 
(/o,0), the distribution predicted by Lynden-Bell is sim- 
ilar to the Fermi-Dirac distribution in quantum mechan- 
ics. When applied to the HMF model, the Lynden-Bell 
theory predicts an out-of-cquilibrium phase transition be- 



tween a magnetized and a non-magnetized phase, and a 
phenomenon of phase re-entrance in the (/o, E) plane [18'. 
The nature of these phase transitions has been studied in 
detail in |19|20|21|22j . In particular, there exist a tricrit- 
ical point separating first and second order phase transi- 
tions, and a critical point (associated with the re-entrant 
phase) marking the onset of a second order azeotropy 
|22) . Direct numerical simulations ^21.231241251261 showed 
a good agreement with the Lynden-Bell prediction in cer- 
tain case^but also evidenced discrepancies in other cases. 
For example, in |21| . the re-entrant phase predicted from 
the Lynden-Bell theory in a very small range of parame- 
ters is confirmed (which is a success of the theory), but 
a secondary re-entrant phase that is not predicted by the 
Lynden-Bell theory is also observed (this secondary re- 
entrant phase has been recently confirmed by another 
group [27] suggesting that it is not a numerical arti- 
fact). More generally, the adequacy, or inadequacy, of the 
Lynden-Bell theory to predict the magnetization of the 
QSS can be read from the numerical phase diagrams plot- 
ted in [HOg. 

These discrepancies can be interpreted as a result of in- 
complete relaxation 15 18 28 29 30". Indeed, the Lynden- 
Bell statistical theory is based on an assumption of ergod- 
icity or, at least, efficient mixing. If the system does not 
mix well, the Lynden-Bell prediction fails and the system 
may be trapped in a steady state of the Vlasov equation 
that is not the most mixed state. This is precisely what 
happens in the case Mq = 1 considered by Antoni & Ruffo 
[8| and Latora et al. |ll|12j . As discussed in |18|29j . the 
failure of the Lynden-Bell prediction is particularly clear 
in that case. Indeed, for an initial condition with Mq — 1, 
we are in the non-degenerate (dilute) limit of the Lynden- 
Bell theory. In this limit, the Lynden-Bell distribution re- 
duces to the Boltzmann distribution (with a different in- 
terpretation). Therefore, as argued in the theoretical 
Boltzmann caloric curve plotted by Antoni & Ruffo [5] and 
Latora et al. |ll|12j should be interpreted as the theoreti- 
cal Lynden-Bell caloric curve. Consequently, the observed 
discrepancies between this theoretical caloric curve and 
the numerical results reveal the failure of the Lynden-Bell 
prediction in that case (this is confirmed by the phase dia- 
grams of |21l26j ). Therefore, the Lynden-Bell theory does 
not explain everything. The limitations of the Lynden-Bell 
theory were emphasized in |18l28j . 

Following these observations, it has been proposed 
[10118131] that Tsallis (/-distributions may provide a good 
fit of the QSSs in certain cases of incomplete relaxatiorj^ 
At the same time, it has been emphasized that this good 
agreement is not expected to be general, i.e. the Tsal- 



We shall present in Fig. 36 a new simulation showing a 



^ This model was first introduced in 1982 by Messer & Spohn 
|3| who called it the cosine model. It was re-introduced in the 
1990s by several authors |4l5l6l7li?] and considerably studied 
since then (see |9I10| for a short historic of the HMF model). 

^ In 2D turbulence, this is called the Miller-Robert-Sommeria 
(MRS) theory [l6lT7] . 



perfect agreement with the Lynden-Bell theory. 

* This is similar to the original claim of Latora et al. [11112) . 
except that we consider incomplete relaxation towards the 
Lynden-Bell distribution (coUisionless regime) , not towards the 
Boltzmann distribution (collisional regime). In the first case, 
the mixing is due to mean field effects, while in the second 
case it is due to discreteness (finite A'^) effects. This is physi- 
cally very different [28j . 
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lis distributions are not universal attractors. Actually, q- 
distributions correspond to what have been called stel- 
lar polytropes in astrophysics ,32j . They were introduced 
long ago by Eddington [33] as particular stationary solu- 
tions of the Vlasov equation. They were used to construct 
simple self-consistent mathematical models of galaxies. 
At some time, they were found to provide a reasonable 
fit of some observed star clusters, the so-called Plummer 
[34] model. Improved observation of globular clusters and 
galaxies showed that the fit is not perfect and more realis- 
tic models have been introduced since then [32 . However, 
stellar polytropes are still important in astrophysics for 
historical reasons and for their mathematical simplicity. 
Similarly, we believe that these distributions will play a 
useful role in the HMF model. Of course, the relevance 
(or irrelevance) of ^-distributions can only be assessed by 
the results of numerical simulations that we now briefly 
review. 

Campa et al. [35j performed direct numerical simula- 
tions of the HMF model starting from an initial condi- 
tion with magnetization Mq = 1 and a waterbag velocity 
distribution corresponding to an energy E = 0.69. They 
obtained a non-magnetized QSS with a velocity distribu- 
tion that they called "semi-ellipse'j^ Chavanis (29j noted 
that this distribution is a particular polytropic (Tsallis) 
distribution with index n — 1. This was the first clear 
evidence of a polytropic QSS in the HMF model. Further- 
more, this distribution has a compact support which is 
very natural in the phenomenology of incomplete violent 
relaxation. Chavanis & Campa [TD] developed a general 
theory of polytropic distributions in the context of the 
HMF model. In this approach, polytropic distributions 
are interpreted as particular steady states of the Vlasov 
equation, like in astrophysics [31]. They studied the dy- 
namical stability problem by using a "thermodynamical 
analogy" and evidenced a rich variety of phase transitions 
depending on the polytropic index n (or q). They com- 
puted the physical caloric curves T^-iniE) and found that, 
for 0.54 < n < 3.56, these curves display a region of nega- 
tive kinetic specific heat Ckin < (see in particular Figure 
23 of [ID]) - They proposed that these results could help in- 
terpreting the "anomalies" reported by Antoni & Ruffo [H] 
that have never been explained so far. 

Pakter & Levin [36] performed direct numerical sim- 
ulations of the HMF model starting from a rectangular 
waterbag distribution with Mq = 0.40 and different val- 
ues of the energy. They found that the QSS has a core- 
halo structur^ The core corresponds to a completely de- 



^ This distribution function differs from the one obtained 
by Latora et al. [12] for the same initial condition. However, 
Campa et al. [35] showed that the ordinary waterbag initial 
condition leads to the presence of large sample to sample fluc- 
tuations so that many averages are necessary. They argued 
that Latora et al. [12] may not have used sufficient averages, 
and they proposed to use isotropic waterbag distributions to 
reduce the fluctuations. 

® This core-halo structure has been also observed in early 
numerical simulations of the Vlasov equation in ID and 2D 
gravity [37138139140141142143144145146] . 



generate distribution in the Lynden-Bell theory (i.e. the 
ground state at T = 0). The halo is interpreted in terms 
of a parametric resonance. This core-halo structure is not 
consistent with the Lynden-Bell prediction that leads to a 
partially degenerate distribution without core-halo struc- 
ture. In addition, for Mq = 0.40, the Lynden-Bell theory 
predicts a second order phase transition T9J while Pak- 
ter & Levin 36J find a first order phase transition. We 
note that the homogeneous core, interpreted as a com- 
pletely degenerate Lynden-Bell distribution, is a polytrope 
n = 1/2 (waterbag distribution) 47 . Therefore, the first 
order phase transition reported by Pakter & Levin |36j 
may be connected to the first order phase transition found 
by Chavanis [U] for the pure waterbag distribution (com- 
pare Figure 6 of [47] to Figure 2 of [^|^ 

Morita & Kaneko [48j performed direct numerical sim- 
ulations of the HMF model starting from initial conditions 
in which the angles and the velocities of the particles have 
Boltzmannian distributions with different temperatures. 
They found initial conditions for which the system does 
not relax toward a QSS. In their simulations, the mag- 
netization exhibits persistent oscillations whose duration 
diverges with N. This long-lasting periodic or quasi peri- 
odic collective motion appears through Hopf bifurcation 
and is due to the presence of clumps (high density regions) 
in phase space. 

This brief review of numerical results in the HMF 
model shows that the nature of the QSS crucially de- 
pends on the initial condition. The purpose of this pa- 
per is to investigate other initial conditions that have not 
been studied previously. We first consider non-magnetized 
initial states with a velocity distribution that is (i) gaus- 
sian, (ii) semi-elliptical, and (iii) waterbag. In each case, 
we determine the physical caloric curve Ti^in (E) of the cor- 
responding QSS and plot the distribution function f{0,v) 
as a function of the individual energy e — + 'P{9). A 
steady state of the Vlasov equation is characterized by a 
function / = /(e). In each case, we find that the caloric 
curve presents a region of negative kinetic specific heat: 
Ckin = dE/dTkin < 0. In that region, we find that the 
QSS is described by a distribution function / = /(e) that 
can be conveniently fitted by a pure polytrope with index 
(i) n = 2, (ii) n = 1, and (iii) n — 1/2. Using the values 
of these indices, we are able to determine the theoretical 
caloric curve T^m (E) and explain the negative kinetic spe- 
cific heat region found in the numerical simulations. For 
lower energies, we find that the system takes a "core-halo" 
structure. The core corresponds to the pure polytrope dis- 
cussed above but it is now surrounded by a halo of parti- 
cles. In case (iii), the distribution function in the core is 
constant (n = 1/2 polytrope) and we recover the results 
of Pakter & Levin [36^ . We also investigate the time evolu- 
tion of the magnetization in the timescale corresponding 
to the QSS. We find persistent oscillations similar to those 



^ We emphasize, however, that the two results are indepen- 
dent since there is no core-halo state in the study of [47j . We 
only suggest that the nature of the phase transition (first or- 
der) is principally due to the n — 1/2 polytropic component 
(the core in [36]). 
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reported by Morita & Kaneko [48]. In cases (i) and (ii), 
they are slowly damped, and in case (iii) they are more 
persistent. We interpret these oscillations in relation to 
the Landau damping around a spatially inhomogeneous 
QSS. For the spatially homogeneous waterbag distribu- 
tion, we know that there is no Landau damping. By con- 
tinuity, the damping rate should be small in the inhomoge- 
neous case, and this may explain the numerical results. We 
also consider unsteady initial conditions with Mq — 1 and 
isotropic waterbag velocity distribution. In that case, we 
find that the system forms phase space holes and dense 
filaments that persist for a very long time. Ultimately, 
these holes disappear and a QSS is formed. We show that 
the kinetic caloric curve is approximately constant]^ cor- 
responding to a polytrope with index no — 3.56. This 
sensibly differs from the results reported in |8lllj . Finally, 
we consider the collisional evolution of an initially Vlasov 
stable distribution, and show that the time-evolving dis- 
tribution function f{v,t) can be fitted by a sequence of 
polytropic distributions with a time-dependent index n{t) 
both in the non-magnetized and magnetized regimes. 

Before discussing these numerical results (Sections ^ 
[6|, we summarize our theory of polytropes jlOI . and pro- 
vide additional theoretical results that were not given in 
our previous paper (Section|3|. In particular, we analyze in 
more detail the nature of the kinetic caloric curve Tkm {E) 
as a function of the polytropic index n since these curves 
are of considerable importance to interpret the numerical 
simulations. Some readers may skip this theoretical part 
and go directly to Section |4] where the results of numerical 
simulations are presented. 



The energy per particle E = H/N can be written as 

i=l 

where M = (Mj -I- Af^)^/^ is the modulus of the magne- 
tization. The equations of motion ([T]) can be written in 
terms of the magnetization components and Aly as 



= -Mj; sin 01 + My cos 0i . (6) 

Following the Kac prescription [49_, the interaction 
strength has been rescaled by a factor 1/N. This is the 
right scaling to properly define the thermodynamic limit 
N — > +00 of a system with long-range interactions such 
as the HMF model. In this way, the total energy of the 
system E ^ N is extensive. However, the energy remains 
fundamentally non-additive [2]. 

The HMF model can be thought of as a set of TV glob- 
ally coupled rotators or XF-spins {6i represents the ori- 
entation of the i-th rotator and Vi is the conjugated mo- 
mentum). In this respect, it is similar to the XY model 
though the interaction is extended to all couples of spins 
instead of being restricted to nearest neighbors. It is also 
similar to a one dimensional periodic self-gravitating sys- 
tem where the potential of interaction has been truncated 
to the first Fourier mode. 



2 The HMF model 

2.1 The Hamiltonian equations 

The HMF model can be viewed as a collection of N parti- 
cles of unit mass moving on a circle and interacting via a 
cosine binary potential 2 . The dynamics of these particles 
is governed by the Hamilton equations 

de, _ dH dvi _ dH 

N N 

where 9i e [0, 27r[ is the angle that particle i makes with 
an axis of reference Ox and Vi = dOi/dt €] — oo,-|-cx)[ is 
its velocity. The HMF model conserves the energy H and 
the number of particles N . The order parameter is the 
magnetization M = {Mx,My) with components 

i i 

* Actually, this curve presents an unexpected hump with a 
region of negative kinetic specific heat. 



2.2 The Liouville equation 

We introduce the A^-body distribution function 
PNiOijVi, ...,9N,VN,t) giving the probability density 
of finding, at time t, the first particle with position 6i and 
velocity Vi , the second particle with position 02 and veloc- 
ity V2 etc. It is normalized such that J Pn Yii dOidvi — 1. 
For an isolated Hamiltonian system, such as the HMF 
model, the evolution of the A'^-body distribution function 
is governed by the Liouville equation 




where Fi — ~ ~W sin(0i — 9j) is the force acting 
on particle i due to the interaction with the other parti- 
cles. It can be written Fi — —MxSm6i + MyCosdi. The 
Liouville equation ([t]) contains exactly the same informa- 
tion as the A^-body Hamiltonian system ([T]). 

2.3 Collisionless evolution: The Vlasov equation 

For systems with long-range interactions, it can be shown 
that the mean field approximation is exact at the ther- 
modynamic limit A^ —J- +00 [50]. This means that the 
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A^-body distribution function is a product of N one-body 
distributions 

PN{dl,Vi, ...,6n,Vn, t)=l[Pl{0^,v,,t). (8) 



Let us introduce the one-body distribution function 

f{o.v,t) = ±{j2fL^sie - e,)S{v - V,)) = p,{e,v,t). it 

is normalized such that 



/[/] 



fdedv = 1. 



(9) 



Substituting the decomposition ([8| in the Liouville equa- 
tion Q and integrating over all variables except ^i, ui, 
we find that, for a fixed interval of time [0,T] (any) 
and N — >■ +oo, the evolution of the distribution function 
f{9,v,t) is governed by the Vlasov equation 



where 



5/ dl d^dl^ 
dt 00 do dv 



<p{e,t)^ / [i-cos{9-e')]p{9\t)de', 

Jo 



(10) 



(11) 



is the self-consistent potential generated by the density of 
particles p{0,t) = J f{d,v,t)dv. The density of particles 
is defined by p{e,t) = ^(E^i ^{0 - Oi)) = Pi{e,t) and 
it is normalized such that J pd0 — 1. The mean force 
acting on a particle located at 6 is {F){9,t) — —d<l>/d9 — 

- /o sm{e-e^p{e', t) de'. Expanding the cosine function 
in equation (11), we obtain 



2.4 The mean field energy 

In the mean field approximation, the energy of the system 
can be written as 

where E^in is the kinetic energy and W the potential en- 
ergy. Using equations (12)-(14), the potential energy can 



be expressed in terms of the magnetization as 



The average kinetic temperature T]^in is defined by 

1, 
2 

Therefore, the energy can be rewritten 
1^ 1-1/2 



(16) 



(17) 



(18) 



For a fixed value of the energy, the kinetic temperature 
directly determines the magnetization, and vice versa. The 
local pressure is defined by 



<^{e,t) = 1 - M^{t) cos 6* - My{t) sine, 



(12) 



p{9,t) = j fv^dv. 
The kinetic energy can therefore be written 

Ekin = ^ J pd9. 
This expression will be useful in the following. 



(19) 



(20) 



where 

M^{t)= p{e,t) cose de, (i3) 



My{t) = J p{e,t) sine de, (14) 

are the components of the mean magnetization M(t). The 
mean force acting on a particle can be written F(e, t) = 
-M^ (t) sin e + My {t) cos 6*. 

The Vlasov equation is a purely mean field equation 
ignoring collisions (more properly, correlations, graininess, 
finite A^-effects...) between particles. It governs the coUi- 
sionless evolution of the HMF model. 

The Vlasov equation admits an infinite number of sta- 
tionary solutions. Any spatially homogeneous distribution 
function / = f{v) is a steady state of the Vlasov equation. 
On the other hand, spatially inhomogeneous distributions 
of the form / = /(e), where e — v'^ /2 + 'P{e) is the individ- 
ual energy, are also steady states of the Vlasov equation. 



2.5 Incomplete violent relaxation 

Starting from an unstable or unsteady initial condition, 
the Vlasov equation is expected to reach, on a coarse- 
grained scale, a QSS. This is called weak convergence in 
mathematics. This QSS is a stable steady state of the 
Vlasov equation. Since the Vlasov equation admits an in- 
finity of steady states, the prediction of the QSS actually 
reached by the system is difficult. A prediction can be 
made based on Lynden-Bell's statistical theory of violent 
relaxation. However, this theory assumes that the evolu- 
tion of the system is ergodic, or at least that the mixing 
is efficient. There are cases where the Lynden-Bell theory 
gives a good prediction. However, there also exist cases 
where this prediction fails. In case of incomplete relax- 
ation, the system can be stuck in a stable steady state of 
the Vlasov equation that is not the most mixed state, i.e. 
that differs from the Lynden-Bell distribution. In this pa- 
per, we consider a particular class of steady states of the 
Vlasov equation, called polytropic distributions, that may 
arise in case of incomplete relaxation. They are associated 
with the Tsallis "entropy" (in the sense given below). 



6 



A. Campa and P.H. Chavanis: Caloric curves fitted by polytropic distributions in the HMF model 



3 Reminders and complements in the theory 
of polytropes 

The theory of polytropes for the HMF model has been 
developed in |10j . Since this theory is rather rich (sev- 
eral cases arise depending on the polytropic index) and 
not well-known, we provide here a summary of the main 
results. In complement to our previous paper, we (i) use 
more conventional notations, (ii) simplify some important 
formulae, (iii) treat a case that was forgotten, (iv) de- 
scribe in detail the physical caloric curve Tkin{E) that 
will be needed to interpret our numerical results, and (v) 
specifically consider the waterbag distribution which is a 
particular polytrope of index n = 1/2. 



3.1 Polytropic distributions in phase space 



The Tsallis entropy is defined by 
1 



S[f] 



iP - f) dedv, 



(21) 



with q > 0. For q — > 1, it reduces to the Boltzmann en- 
tropy 

S[f] = - f flnfdddv. (22) 



We consider the microcanonical problem 



max {S[f] 



E[f]^E, I[f] = l}, 



and the canonical problem 



mm 
/ 



{F[f]=E[f]-TS[f] I /[/] = 1}. 



(23) 



(24) 



As explained in |10l51j . these optimization problems de- 
termine particular steady states of the Vlasov equation of 
the form / = /(e) with /'(e) < that are dynamically 
stable. In this dynamical interpretation, S* is a pseudo en- 
tropy, f is a pseudo free energy, and T > is a pseudo 
thermodynamical temperature. By an abuse of language, 
and to simplify the terminology, we shall omit the prefix 
"pseudo" in the sequel. However, we must keep in mind 
that all the references to thermodynamics in our study are 
purely effective: we are dealing with dynamical stability, 
not thermodynamical stability. Nevertheless, it is conve- 
nient to develop a thermodynamical analogy and use a 
similar terminologjj^ 



We recall that canonical stability (criterion (24|) im- 
plies microcanonical stability (criterion (pSl)) but the con- 
verse is wrong in case of ensembles inequivalence i53i . Fur- 
thermore, we recall that the thermodynamic-looking op- 
timization problems (23) and (24) provide just sufficient 



conditions of Vlasov dynamical stability. More refined dy- 
namical stability criteria are discussed in |51j . 



® Tsallis and co-workers have tried to give a real thermo- 
dynamical interpretation to the functional (21 1. Here, we just 



consider its Vlasov dynamical stability interpretation [10151] 
and refer to [31152] for other interpretations. 



The critical points of entropy at fixed energy and nor- 
malization are determined by the variational principle 



SS ~ (3SE ~ aSI = 0, 



(25) 



where /3 = 1/T and a are Lagrange multipliers (T is the 
thermodynamical temperature and a the chemical poten- 
tial). This yields the Tsallis g-distributions 



f(0,v) 



l/(-3-l) 



(26) 



where fJ. ^ [1 — (g — l)a\/q. The notation [x]^ stands 
for [a;]+ = .T if X > and [x]-^- = if a; < 0. These 
distributions are also obtained as critical points of the free 
energy at fixed normalization, satisfying SF+aTSI = 0. In 
astrophysics, they are known as "polytropic distributions" 
[32j . They are generally labeled by a polytropic index n 
that is related to the parameter q by the relation 



1 



(27) 



For n = 1/2 (q — > +oo), the polytropic distribution func- 



tion (26 1 reduces to the waterbag, or Fermi, distribution 
(see Section 3.8). For n +oo {q — >• 1), it reduces to the 



isothermal (Maxwell-Boltzmann) distribution 



f{e,v)^J^e 



-/3[i#+*(e)] 



(28) 



We need to distinguish two cases depending on the sign 
of q — 1. (i) For q > 1 (n > 1/2), the distribution function 
can be written as 



/ ^ A{era - e)] 



(29) 



where we have set A — [f5{q — l)/q]'^ and e,„ = 
q^i/[l3{q — 1)]. Such distributions have a compact sup- 
port in phase space since they vanish for e > em- At 
a given position 6, the distribution function vanishes at 
V = v,n{9) = y/2{e,n-<P{e)). For n = 1/2, e„ and vUO) 
correspond to the Fermi energy and to the Fermi velocity, 
respectively, (ii) For < g < 1 (ri < —1/2), the distribu- 
tion function can be written as 



/ 



A 



(eo + e)i-<! 



(30) 



where we have set A = [f5(l — q)/q]^ and eo = qfi/[/3{l — 
q)]. Such distributions are defined for all individual ener- 
gies. At a given position 9, the distribution function be- 
haves, for large velocities, as / ~ y-'^/i'^-i) ^ 
In the following, we shall only consider distribution func- 
tions for which the density and the pressure are finite. 
This implies 1/3 < g < 1 (n < -1). 

Polytropic distributions may arise as a result of incom- 
plete relaxation. In that case, the system does not mix well 
enough to justify the establishment of the Lynden-Bell 
distribution which corresponds to the most mixed state. 
In general, incomplete relaxation manifests itself by the 
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fact that mixing takes place only in a subdomain of phase 
space. This generally leads to distribution functions with a 
compact support (the distribution function vanishes above 
a certain energy em) while the distribution function pre- 
dicted by Lynden-Bell remains strictly positive for all val- 
ues of the energy (which is rather unphysical) . Polytropes 
with index n > 1/2 have this property of confinement 
and this is why, in our opinion, th ey p lay an important 
role in case of incomplete relaxatiorp^ Therefore, in the 
following, we shall essentially consider the case n> 1/2. 



3.2 Polytropic equation of state 

To any distribution function of the form / — /(e) with 
/'(e) < 0, we can associate a corresponding barotropic gas 
by defining the density p = J f dv = p{'P) and the pressure 
P — I /"^^^ — p(^)j and eliminating 'P{6) between these 
two expressions. This determines a barotropic equation of 
state p = p{p)- For the polytropic distribution (26 1, we 

m-. 



get 



with 



p = Kp'^, 7 = H--, 



^^^|y2^r(l/2)r(l/2 + n) 



71 + 1 

for n > 1/2 and 



r{n + i) 



K 



n + l\ r(l/2-n) J 



(31) 



(32) 



(33) 



for n < —1. Equation (31) is the well-known polytropic 
equation of state. This is the reason why the distributions 
(26) are called polytropic distributions. The polytropic 
constant K is sometimes called the "polytropic temper- 
ature". We can show that K \s a. monotonically increas- 
ing function of the thcrmodynamical temperature T 10 . 
For isothermal systems [q — \, n ~ oo, 7 = 1), we have 
K = T. 



3.3 Polytropic distributions in physical space 



The density is obtained by integrating Eq. (26) over the 
velocity. We find that the density is related to tEe potential 
^'(61) by |IOj: 



where A = em/{K{n+l)) for n > 1/2 and A — eq/ {—K{n+ 
1)) for n < — 1. For 7 = 1, Eq. (34) reduces to the isother- 
mal (Boltzmann) distribution 



p(0) = Ae-^*W, 



(35) 



with A = (27r//3)^/^A'. The polytropic distribution in 
physical space p = p{'l>) given by Eq. ( [34| has the same 
mathematical form as the polytropic distribution in phase 
space / = /(e) given by Eq. (26) with 7 playing the role 



of q and K playing the role of T — 1/ /3. In this correspon- 
dence, 7 is related to q by Eqs. (pTI) and (31) leading to 



7 = (3g - l)/(g 1) and K is related to T by Eqs. ([32) 



and ( 33 ) . Polytropic distributions (including the isother- 
mal distribution) are apparently the only distributions for 
which /(e) and p{'P) have the same mathematical form. 
Using Eq. ([34|, the polytropic distribution function 



(26 1 can be rewritten [TU] : 



n-1/2 



(36) 



where Z is given for n > 1/2 by 



and for rt < — 1 by 

For n — >■ -|-oo, we recover the isothermal distribution 

1/2 



(38) 



p{9)e- 



(39) 



Other useful expressions of the polytropic distribution 
function are given in [TU] . 

3.4 Homogeneous phase: Jeans-type instability 

It is convenient to define the polytropic temperature by 



= K 



7-1 



(40) 



P{0) = 



(34) 



+ 



Of course, polytropic distributions are not the only distri- 
butions with a compact support and, indeed, we will find that 
the QSSs may be described by more complicated distributions 
(e.g. core-halo states). In other words, Tsallis distributions are 
not universal attractors in case of incomplete relaxation. How- 
ever, we have suggested |18I10| that polytropic distributions 
may provide a good fit of QSSs in certain cases, and we will 
give numerical evidence of that claim in Sections [4]|6| 



In the homogeneous phase (M = 0), using Eqs. (36) and 
( 40 1 , the polytropic (Tsallis) distributions can be written 



fiv) = C 



1 - 



2{n + l)e 
where C is given for n > 1/2 by 



-1/2 



(41) 



1 r(n+l) 1 
27r r{l/2)r{l/2 + n) ^2{n + 1)0 ' ^ ^ 
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and for n < — 1 by 

1 r(l/2-n 



1 



27rr(l/2)r(-n) ^-2(n + l)©' ^ ^ 

For n > 1/2, the maximum v eloc ity is Vmax — 
i/2(ri + 1)6*. For n +oo, Eq. (41) reduces to the 
isothermal (Maxwell) distribution 



/(«) = 



1/2 



pe 



(44) 



On the other hand, using Eqs. (201, (31 1 and (40), we find 



that the energy can be written as 



(45) 



Therefore, in the homogeneous phase, the kinetic temper- 
ature coincides with the polytropic temperature: 



Thin — O. 



(46) 



In addition, the caloric curve 0{E) is the same for all 
indices n. It is defined for Tkin = > and E > 1/2. 
The specific heat is C = dE/dO = dE/dTkin = 1/2. The 
equality ( 46 1 is not true anymore in the inhomogeneous 
phase. 



For the polytropic equation of state (31 1, the square of 



the velocity of sound in the homogeneous phase is given 

by 



70. 



(47) 



It can be shown that a spatially homogeneous distribu- 
tion function f{v) with a single maximum at t; = is 
dynamically stable if, and only, if 



1 



-|-CX3 



fiv] 



dv > 0. 



(48) 



For a distribution function of the form / = /(w^) with 
/'(w^) < 0, the stability criterion can also be written as 



(49) 



These two stability criteria are equivalent and they de- 
termine the spectral, linear and formal stability of the 
distribution [,9 51.54.55) . According to these stability cri- 
teria, the homogeneous phase of a polytropic distribution 
is stable if [5110] : 



e>er = 



27' 



E> Er = 



1 

47 



1 

2' 



(50) 



and unstable otherwise. For isothermal distributions (n = 
oo,7 = l,q = l),we recover the critical values (Tc, Ec) ~ 
(1/2, 3/4). For the waterbag distribution (n = 1/2, 7 = 3, 
q = +00), we get (6>c,£'c) = (1/6,7/12). For the semi- 
ellipse (n = 1, 7 = 2, g = 3), we obtain {0c, Ec) = 
(1/4,5/8). As shown in previous works [51619155] . the in- 
stability below Ec is similar to the Jeans instability in 
astrophysics [32j . 



3.5 Inhomogeneous phase: complete and incomplete 
polytropes 

We now consider spatially inhomogeneous polytropic dis- 
tributions. We can assume without loss of generality that 
the distribution is symmetrical with respect to the x- 
axis {i.e. with respect to the angle f? = 0). In that case, 
the potential can be written (l>{0) = 1 — Mcos9 where 

M = p{9) cos 9 d9 is the magnetization {My = and 
M = Mx). The density profile (34) can be rewritten [TU] : 



p{e) = A[k + 



cos 9 



where we have noted 



M 



(51) 



(52) 



We must consider two cases k = ±1. As realized in [47] , 
the case k = —1 was forgotten in [10]. However, this for- 
getfulness does not alter the study of phase transitions 
made in [TO] since the branch k = — 1 only completes the 
caloric curve up to the ground state E ~ 0. 

For X > 0, the density profile is concentrated around 
^ = and for a; < 0, we get a symmetrical density pro- 
file concentrated around 9 = n. We can therefore re- 
strict ourselves to a; > 0. For x = 0, we recover the 
homogeneous distribution p — l/(27r). For x > 0, the 
density profile is monotonically decreasing. Let us first 
assume n > 1/2 (i.e. 1 < 7 < 3) and k = -1-1. If 
X < Xc = n + 1 = 7/(7 — 1), the density is strictly positive 
on — TT < 9 < TT (incomplete polytrope). If x > Xc, the 
density has a compact support (complete polytrope). In 
that case, it vanishes for \9\ > 9c ~ arccos(— Xc/x). We 
note that 9c > tt/2. For x +00, the density profile is 
given by ^(6*) = r{{2 + n)/2)/[y^r{{l + n)/2)]cos'' 9 
We now assume n > 1/2 and k — —1. The central density 
is defined only for x > Xc. In that case, the polytrope is 
complete: the density vanishes for \9\ > 9c = arccos(xc/x). 
We note that 9c < tt/2. For x — t' Xc, the distribution tends 
to a Dirac peak p{9) — 6{9). Finally, we consider the case 
n < —1 (i.e. < 7 < 1) and k — +1. The central density 
is defined only for x < Xc = |n -|- 1| = 17/(7 — 1)| and the 
polytrope is incomplete. Some typical density profiles are 
represented below in Figures [18] and [19] 



The amplitude A of the density profile (51) is deter- 
mined by the normalization condition J pd9 — 1, leading 
to 

A= ^ / . ^ . (53) 
27r/^.o(x) ^ ' 

where we have introduced the 7-deformed modified Bessel 
functions ^0^: 



(x) 



1 

2^ 



7- 



7 



X cos 9 ) cos{m9) d9. 



^ (54) 
On the other hand, substituting Eq. (51) in Eq. (13) for 
the magnetization, we obtain the self-consistency relation 



M 



ly^l{x) 



(55) 
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Combining Eqs. ( 52 1 , ( 53 1 and ( 55 1 , wc find that the poly- 
tropic temperature ( 40 ) is given by 



^ 1 -^7,1 (^) 

a;/^,o(a;)2-T- 



(56) 



Finally, after some calculations [TD] , we can show that the 
kinetic temperature and the energy are given by 



M 



Tkm = - — + n— 
7 X 



E 



1 - M2 



(57) 



(58) 



Equation (57) is equivalent to equation (109) of (TUl af- 
ter a straightforward simplification. From equations (551- 
(581), we can obtain the curves 0{E), M{0), M{E), and 



Tkin{E) in parametric form with parameter x. Some of 
these curves have been analyzed in detail in [10 . This 
analysis will be completed in Section [3. 6| 

Close to the bifurcation point (i?c, Oc)^ i-e. for a; — > 0, 
using the asymptotic expansions given in [10 , we find that 
the thermodynamical specific heat C = dE/dO is given by 



C 



- 57 - 2 



2(2-7) 

Let us introduce the critical indices 



7* 



41 



1 



2.8507811. 



0.54031242. 



(59) 



(60) 



(61) 



The specific heat of the inhomogeneous phase at the bifur- 
cation point is positive for n > 1, negative for < n < 1 
and positive again for 1/2 < n < . It is infinite for n — 1 
and vanishes for n — n^,. For n < — 1, the specific heat is 
always positive. These results explain the behavior of the 
caloric curve 0{E) close to the bifurcation point ^OJ. For 
the waterbag distribution {n = 1/2, 7 = 3, g = -t-oo), we 
get C = 1/2 like in the homogeneous phase [47]. For the 
semi-ellipse (n = 1, 7 = 2, g = 3), we obtain C = 00. On 
the other hand, close to the critical point, the magnetiza- 
tion is given as a function of the energy by 



87 



272 - 57 



(62) 



For isothermal distribution {n 



-00), the density 



profile (51) takes the form 



p{0) 



1 



2TrIo{x) 



^PM cos 8 



(63) 



with X — M/T. The self-consistency relation ( 55 ) reduces 
to 

(64) 



Finally, the temperature ( 56 ) or ( 57 1 and the energy ( 58 ) 
become 



T — Tkin — 



M 



1 1 - A'P 



(65) 



(66) 



Close to the bifurcation point {Ec,Tc)={l/2, 3/4), the 
specific heat is C = dE/dT — 5/2 and the magnetization 
is M = y/8/5{Ec - £^)^/2 Tiiis returns the well-known 
results of the Boltzmann thermodynamical analysis (see, 
e.g., i). 



3.6 The physical caloric curve 

In our previous paper |10| , we have plotted the thermody- 
namical caloric curves giving the thermodynamical tem- 
perature T, or the polytropic temperature (9, as a func- 
tion of the energy E. This is the correct way to determine 
the stability of the system in canonical and microcanoni- 
cal ensembles and study phase transition^^ However, the 
temperature that is directly accessible to the experiments 
or to the numerical simulations is the kinetic temperature 
Tkin- In general, T^m 7^ T = 1//3 in a QSS. For compari- 
son with the numerical results of Sections |4]|6j it is useful 
to study the physical caloric curve Tf^iniE) as a function 
of the polytropic index n. This study has been only partly 
done in our previous paper (see Figures 6 and 23 of |10j ) 
and it is here systematically continued. 

Close to the bifurcation point {Ec, T^in): i.e. for x ^ 0, 
using the asymptotic expansions given in 10^ , we find that 
the physical (or kinetic) specific heat Ckin = dE/dTkin is 
given by 



272 - 57 - 2 



2(272-7-2)- 
Let us introduce the critical indices 



70 



no 



1 + 717 



17-3 



1.2807764.. 



3.5615528. 



(67) 



(68) 



(69) 



The kinetic specific heat of the inhomogeneous phase at 
the bifurcation point is positive for n > no, negative for 
n* < n < no and positive again for 1/2 < n < n*. It is 
infinite for n = tiq and vanishes for n — n^,. For n < — 1, 
the kinetic specific heat is always positive. These results 
explain the behavior of the physical caloric curve Tkin[E) 
close to the bifurcation point (see Figure [l]). We note that 
for 1 < n < no, the solutions close to the bifurcation point 



We recall again that we are actually considering the Vlasov 
dynamical stability of polytropic distributions using a thermo- 
dynamical analogy. 
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Fig. 1. The physical caloric curve for different values of n. 




Fig. 2. The physical caloric curve for n > no (specifically 
n = oo). The kinetic specific heat is positive. The physical 
caloric curve exhibits a second order phase transition at Ec- 




Fig. 3. The physical caloric curve for n = no — 3.56. The 
kinetic specific heat close to the bifurcation point is infinite. 



have negative kinetic specific heat Ckm = dE/dTkin < 
although they have positive thermodynamical specific 
heat C = dE/dT > \10\. Therefore, one may observe 
negative kinetic specific heats although the ensembles are 
equivalent. A similar observation has been made in |56j 
for the Lynden-Bell distribution. 



0.8- 



0.6 



0.4- 



0.2- 



E' T, . 

ki 

u 



0.2 



0.4 



0.6 



0.8 



Fig. 4. The physical caloric curve for umcp — 0.68 < n < no 
(specifically n = 1.5). The kinetic specific heat close to the 
bifurcation point is negative although the thermodynamical 
specific heat is positive |10| . 




Fig. 5. The physical caloric curve for n = 1. The kinetic spe- 
cific heat close to the bifurcation point is negative (Ckin ~ 
— 1/2). For 9/16 < E < Ec = 5/8, the physical caloric curve is 
a straight line given by T^in = 3/2 — 2i5. The thermodynamical 
specific heat is infinite |10) . 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Fig. 6. The magnetization curve for n > umcp (specifically 
n = 1). For 9/16 < E < Ec = 5/8, we have M = 2{Ec-Ef/'^. 
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Fig. 7. Microcanonical phase diagram for n > umcp — 0.68. 
It shows a second order phase transition at and a region of 
negative kinetic specific heat between E' and Ec for umcp < 
n < no. 




Fig. 8. Evolution of r^j„ and T^i„ as a function of n. 



We must be careful that we cannot deduce any stability 
result in the canonical ensemble (the optimization prob- 
lem (24)) from the physical caloric curve Tkin{E) since 



Tfcin is not the correct variable (the correct variable in the 
canonical ensemble is T or 10 ). Therefore, we restrict 
ourselves to the microcanonical ensemble (the optimiza- 
tion problem (23l). Actually, this is the most appropriate 
ensemble since the natural control parameter is the energy. 
Since phase transitions for polytropic distributions have 
already been discussed in our previous paper (from the 
thermodynamical caloric curve 0{E)), our discussion here 
will be more concise. We shall denote global entropy max- 
ima by (S), local entropy maxima by (M) and minima or 
saddle points of entropy by (U). In the thermodynamical 
analogy, these symbols correspond to stable, metastablep^ 
and unstable states. Coming back to the dynamical inter- 
pretation, stable states (S) and metastable states (M) are 
dynamically Vlasov stable. We cannot definitely conclude 



For systems with long-range interaction, the metastable 
states are extremely robust. In practice, they should be con- 
sidered on the same footing as stable states. However, in our 
theoretical analysis, we shall distinguish between fully stable 
and metastable states. 




Fig. 9. The physical caloric curve for umtp — 0.563 < n < 
nMCP (specifically n — 0.6). 




Fig. 10. Enlargement of the physical caloric curve for umtp < 
n < Umcp (specifically n = 0.6). The kinetic specific heat 
close to the bifurcation point is negative. The kinetic caloric 
curve displays a second order phase transition at Ec and a first 
order phase transition at Et . The transition energy Et has been 
determined from the entropy curve S{E) as explained in [10) . 
Note that the Maxwell construction cannot be performed on 
the physical caloric curve Tkin{E) since Tkin i^T ~ 



that unstable states (U) are Vlasov unstable, except for 
homogeneous distributions, since the optimization prob- 
lem ( |23[ ) provides just a sufficient condition of dynamical 
stability [5l\. 

For n > Umcp — 0.68, the physical caloric curve 
TkiniE) exhibits a second order phase transition marked 
by the discontinuity of dTkin/dE at E = (see Figures 
[2][3j|4j andjsj. An example of magnetization curve is given 
in Figure 161 For n > uq ~ 3.56, the kinetic specific heat 
is positive (see Figure [2]). For umcp < n < uq, there is 
a region of negative kinetic specific heat between E' and 
Ec (see Figures [4] and [5| . The microcanonical phase di- 
agram in the (n,E) plane for n > umcp is represented 
in Figure 171 For completeness, we have also plotted T^j„ 
and T^j^^^ (the latter being the kinetic temperature corre- 
sponding to E') as a function of n in Figure [sj although 
Tkin should not be regarded as a control parameter. For 
Umtp — 0.563 < n < umcp, there is a very interesting 
situation, already noted in IlOj, in which the caloric curve 
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Fig. 11. The magnetization curve for umtp < n < umcp 
(specifically n = 0.6). 




Fig. 14. The physical caloric curve for n = n* ~ 0.54. The 
kinetic specific heat close to the bifurcation point vanishes. 
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Fig. 12. The physical caloric curve for n, < n < umtp (specif- 
ically n = 0.55). The kinetic specific heat close to the bifurca- 
tion point is negative. There is a first order phase transition at 
E = Et. The metastable branch exhibits a second order phase 
transition at Ec- 



n=0.55 



Fig. 15. The physical caloric curve for 1/2 < n < n* (specifi- 
cally n = 0.5). The kinetic specific heat close to the bifurcation 
point is positive. There is just a first order phase transition. 
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Fig. 13. The magnetization curve for n* < n < umtp (specif- 
ically n = 0.55). 



Fig. 16. The magnetization curve for 1/2 < n < n* (specifi- 
cally n = 0.5). 
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Fig. 17. Microcanonical phase diagram. 



exhibits a second order phase transition at Ef. between 
the homogeneous phase and the inhomogeneous phase (as 
before) and a first order phase transition at E t between 
two inhomogeneous phases (see Figures 9 f 1 ) . The first 
order phase transition is marked by the discontinuity of 
Tkin aX E = Et and the existence of metastable branches. 
Therefore, umcp — 0.68 and Emcp — 0.5873 is a mi- 
crocanonical critical point marking the appearance of the 
first order phase transition. At n = umtp, the energies 
of the first and second order phase transitions coincide 
{Et = Ec) and, for 1/2 < n < umtp, there is only a 
first order phase transition aX E = Et between homoge- 
neous and inhomogeneous states (see Figures [12] [l3j [l4j 
and 



15). Therefore, umtp ^ 0.563 and Emtp ^ 0.5901 is 



a microcanonical tricritical point separating first and sec- 
ond order phase transitions. For ~ 0.54 < n < umtp, 
there remains a sort of second order phase transition at 
E = Ec for the metastable states (see Figures 12 and 13). 
Between E'^ and Ec, the kinetic specific heat is negative. 
At n = n*, the kinetic specific heat close to the critical 
point vanishes and the energies E'^ and Ec coincide (see 
Figure 14). For 1/2 < n < n*, the kinetic specific heat 
close to the critical point is positive and the "metastable" 
second order phase transition disappears (the inhomoge- 
neous states close to the bifurcation point are unstable). 
In that case, there is only a first order phase transition 
(see Figures [Ts] and 16). The microcanonical phase dia- 
gram in the Xji,E) plane summarizing all these results is 
represented in Figure [iTj 

In conclusion, for n > umcp we have second order 
phase transitions, for umtp < n < timcp we have first 
and second order phase transitions, and for 1/2 < n < 
Umtp we have first order phase transitions. 

Remark 1: As a corollary, we emphasize that for n > 
nucp — 0.68, the inhomogeneous polytropes are always 
entropy maxima at fixed energy and normalization so they 
are dynamically Vlasov stable. This is an important the- 
oretical result of ilOi. 



3.7 The polytrope n = l 

The polytropic index n — \ (i.e. 7 = 2, g = 3) is particular 
because it corresponds to a canonical tricritical point [10' . 



Furthermore, for this index, the algebra greatly simplifies 
and analytical results can be obtained. This is because the 
relationship ( 34 ) between the density and the potential is 
lineair3 



The spatially homogeneous distribution is 



1/2 



(70) 



if |w| < Vmax = 2-\/0 and / = otherwise. It corre- 
sponds to what has been called "semi-ellipse" in . The 
homogeneous phase is dynamically stable if, and only, if 
0>ec = ll^OT: E>Ec = 5/8. 

The density profile of n = 1 polytropes is 



p{e) = a(k+|cos6i) 



(71) 



Some density profiles are plotted in Figure [TSj For k = +1 
and X < Xc = 2 (incomplete polytropes), the deformed 
Bessel functions take the simple form l2fi{x) — 1 and 
-^2,1(2^) — x/4,. Then, we get 



A 



E = 



1 

27r' 



M 



X 

4' 



64' ^'=""4 + 32- 



(72) 



(73) 



The spatially inhomogeneous distributions with x < Xc 
have the same polytropic temperature 0c — 1/4 but dif- 
ferent energies ranging from 9/16 to Ec = 5/8. Their 
magnetizations are in the range < AI < 1/2, and they 
are related to the energy by M — 2^/Ec — E (see Figure 
[6|. The thermodynamical caloric curve forms a platearp^ 
0{E) = 0c in the range 9/16 < < 5/8, so it has an infi- 
nite specific heat C = dE/d0 — 00 (see Figure 12 of [10]). 
By contrast, in the same range of energies, the physical 
caloric curve (see Figure™ is given by 



Tkin — 2 ^2-^' 

and it has a constant specific heat 



Cki' 



dE 
dTi 



kin 



(74) 



(75) 



which turns out to be negative. Therefore, for n = 1, the 
thermodynamical and physical caloric curves are very dif- 
ferent. 

For complete polytropes, we can get analytical expres- 
sions of the thermodynamical parameters by using the 



Actually, except for these nice mathematical properties, it 
is not clear whether the polytrope n = 1 plays a special role 
in the physics of the problem. As we shall see in the numerical 
part of the paper, other polytropic distributions are relevant 
as well. This soften the claim made in our previous paper [lUj . 

These solutions are stable in the microcanonical ensem- 
ble and metastable in the canonical ensemble [Tn|. This corre- 
sponds to a situation of partial ensemble equivalence [53] . 
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Fig. 18. Density profile of n = 1 polytropes for different values 
of X. We have added the case k — ~1 that was forgotten in 

noi. 



following expressions of the 7-deformed modified Bessel 
functions: 



^2,0(2;) 



^2,1(2^) 



1 
1 

271 



\/x^ - 4 + 2k 



arccos 



4 + — arccos 



2k 

X 

2k 

X 



(76) 
(77) 



Other analytical results are given in [lOj . 

3.8 The waterbag distribution (polytrope n = 1/2) 

A detailed study of the possibly inhomogeneous waterbag 
distribution that is a steady state of the Vlasov equation 
has been performed in j47i . Here, we only recall the main 
results of this study that will be needed to interpret the 
numerical simulations of Section [4731 The waterbag distri- 
bution is defined by 



(e < ep), 
(e>ef). 



(78) 



It is similar to the Fermi distribution in quantum mechan- 
ics where /o plays the role of the maximum value of the 
distribution function fixed by the Pauli exclusion principle 
and e p plays t he r ole of the Fermi energy. Comparing Eq. 
( 78 ) with Eq. ( 29 ) , we see that the waterbag distribution 
corresponds to a polytrope of index n = 1/2 (i.e. 7 = 3, 
q = -l-oo). This correspondence can also be obtained by 
determining the equation of state associated with the wa- 
terbag distribution. To that purpose, we rewrite Eq. (78 1 
in the form 



(v < vpie)), 

{v>vf{0)), 



(79) 



where vp{6) — yj2{ep — ^(0)) is the local maximum ve- 
locity (the analogous of the Fermi velocity in quantum 
mechanics). The density is given by p{6) = 2fQVp{9) and 



the pressure by p{9) — {2/3)foVp{d). Eliminating vp{6) 
between these two expressions, we find that the equation 
of state of the waterbag distribution is 



1 



P 



12/^' 



(80) 



This is the equation of state of a polytrope of index n = 
1/2 and polytropic constant K = 1/(12/q). Therefore, /o 
is related to the polytropic temperature = K/{2tt)^ by 

fo = ^4^. (81) 



47rV36' 

The velocity of sound is given by 
1 



2/0 



Pi0)^vpie), 



(82) 



and it coincides with the maximum local velocity (the 
Fermi velocity in quantum mechanics). 

The spatially homogeneous waterbag distribution is 
the step function 



/ = /o, if \v\ < vp, 
/ = 0, if \v\>vp, 



(83) 



where vp is determined by the normalization condition ^ 
leading to 



Vp 



47r/o 



The kinetic temperature is given by Tki 
the total energy can be written as 



E = 



(84) 

= vj,/3, so 
(85) 



The velocity of sound in the homogeneous phase is Cs = 
Vp. According to the general criteria (48) and (49), the 
homogeneous phase is dynamically stable if 



E>Ec 



fo < (/o)c = 



7 

' 12' 
1 



0>0c 



1 

6' 



27r%/2' 



Vp > {vp)c 



1 

71' 



(86) 



(87) 



and unstable otherwise. These stability criteria can also 
be directly obtained from the dispersion relation e(aj) = 
which can be solved analytically for the waterbag distri- 
bution (see, e.g., [9]). 

Using the results of Section |3.5[ the density profile 
of the spatially inhomogeneous waterbag distribution is 
given by 



27r/3^o(a;) 



-X cos 6 



1/2 



with K ~ sgn(ei? — 1). The polytrope is incomplete for 
K = +1 and X < Xc ^ 3/2, and complete otherwise. Some 
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Fig. 19. Density profile of n = 1/2 polytropes for different 
values of x. For k = +1 we have taken x — 1 and a:: = 3. For 
K = — 1 we have taken x — 3. 
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Fig. 20. Phase space portrait of n = 1/2 incomplete {x = 1) 
polytropes with n — +1. The distribution function has the 
uniform value / = 770 in the shaded area (e < ep) and f — 
outside (e > ep). 



representative density profiles are plotted in Figure |19| 
The local maximum velocity is related to the local density 
by vp{9) — p{6)/{2fo), so it is proportional to the density 
profile. Combining Eqs. (56), (81), and (88), we obtain 



vf{0) 



'3/3,1(2:) 



2 

-X cos ( 
3 



1/2 



(89) 



, a;/3,o(a;) 

The waterbag distribution is represented in phase space 



in Figures 20 and 21 The velocity distribution (/>(«) can 
be obtained~analytically [17] . 

The dynamical stability of the waterbag distribution 
has been studied in [4T. The magnetization curve M{E) 
is reproduced in Figure [16] and the physical caloric curve 
Tkin[E) is plotted in Figure 15 
transition at Et 



There is a first order phase 
ced by the discontinuity of 
0.44 to zerc{3 The 



0.588 mari^ 
the magnetization passing from Mt 

homogeneous phase is stable for E > Et, metastable for 
Ec = 7/12 < E < Et and unstable for < £; < E^. 
The inhomogeneous phase is stable for {0 < E < Et, 
Mt < M < 1), metastable for {Et < E < E^ = 0.595, 
= 0.37 < M < Mt), and unstable for (E^ < E < E^, 
0<M < M,). 

Remark 2: The magnetization curve M{E) of the pure 
waterbag distribution (polytrope n — 1/2) reported in 
Figure 16 (see also 07]) be compared with Figure 2 
for a core-halo state in which the 
1/2. They both exhibit a first 



of Pakter & Levin [36] 
core is a polytrope n 
order phase transition. 

Remark 3: As discussed in 47 , the waterbag distri- 
bution corresponds to the minimum energy state in the 
Lynden-Bell theory, when the initial condition has only 
two levels and /q. This has been used to construct the 
curve Egrou7id{fo) in Appendix A of [55]. 
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Fig. 21. Phase space portrait of n = 1/2 complete {x — 3) 
polytropes with k — +1. The distribution function has the 
uniform value / = r/o in the shaded area (e < ep) and / = 
outside (e > e^). 



4 Numerical simulations 

We have performed a number of simulations of the HMF 
model that we analyze in this and in the following two 
Sections. We have considered different types of initial con- 
ditions. In this Section, we study the characteristics of the 
QSS that are reached by the system after a violent relax- 
ation from a spatially homogeneous initial state which is 
Vlasov unstable. In Section [5] we study the evolution of 
the system initially prepared in a Vlasov stable spatially 
homogeneous state. In that case, the system is in a QSS 
from the beginning, and the slow evolution is due to finite 
size effects ("collisions"). We have already analyzed this 
evolution in the non-magnetized regime in Ref. |57| . and 
here we extend our analysis to the magnetized regime. In 
Section [6] we treat the case in which the system is initially 

As discussed in detail in [17] 1 the thermodynamical caloric 
curve 0{E) of the waterbag distribution is very particular since 
the temperature jump at the transition point shrinks to zero 
as n — >■ 1/2. It is therefore more convenient to study the phase 
transition on the magnetization curve M{E) or on the physi- 
cal caloric curve TkiniE) than on the thermodynamical caloric 
curve 0{E). 
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in an unsteady inhomogeneous state with magnetization 
Mq = 1 and isotropic waterbag distribution. 

In this Section, the initial distributions fo{v) of the ve- 
locities are of several types: Gaussian, i.e., fo{v) cx e~°" , 
semi-elliptical, i.e., fo{v) ~ v^, and waterbag, i.e., 

/o(w) oc 0{vm — \v\), with the value of the parameters a, 
vq and determining the energy E. These three cases 
correspond to homogeneous polytropic distributions with 
n = oo, n ~ 1 and n = 1/2, respectively. In all cases, 
the energy E is chosen smaller than the critical energy Ec 
for Vlasov stability of the corresponding case. The critical 
energy, computed from Eq. (50), can be written in terms 
of n as 

Substituting the values of n we get that Ec = 0.75 for 
the Gaussian distribution, Ec — 5/8 for the semi-elliptical 
distribution, and Ec = 7/12 for the waterbag distribution. 
The state is unstable if E < Ec- Then, after a rapid relax- 
ation, the system settles down in a QSS. In this Section, 
our purpose is to analyze the physical caloric curves of the 
QSSs and to compare them with those of the polytropic 
distribution functions. 

The simulations have been performed with N — 2^^ 
particles. However, they are representative of a system 
with twice the number of particles, i.e., 2^^, gaining in 
this way a factor of 2 in the computer time needed. In 
fact, we have exploited the following symmetry property 
of the HMF model. If, in the initial conditions, for each 
particle with {9,v) there is a particle with (— 6',— w), this 
property is conserved throughout the dynamics, keeping 
at the same time My = 0. In fact, choosing the numbering 
of the particles in such a way that in the initial conditions 



we have Oi^k = 
motion ^ give: 



-Oj and u, 



the equations of 



-RLr sin 6: 



Mx sin 6., 



that proves the statement. Furthermore, it is sufficient to 
follow (and thus to represent in the computer) the dynam- 
ics of only the first ^ particles. We remark that extracting 
the initial condition from a distribution that is invariant 
under the inversion {9, v) — )> (—0, —v) should not intro- 
duce any peculiarity, since we expect that in the thermo- 
dynamic limit this invariance is satisfied. 

We should note also the following point. The length 
of our simulations is sufficient to follow the rapid relax- 
ation to the QSS and its initial stages. With 2^^ particles 
the lifetime of the QSS is expected to be very large, since 
the collisional processes of the dynamics, due to finite size 
effects, are very slow. It is true that the lifetime of magne- 
tized QSSs, as the ones we are going to analyze, is expected 
to be proportional to N and not to a higher power 
of N as for homogeneous QSSs [35]; however, such high 
values of N lead to lifetimes that are much larger than the 
times we have considered. In conclusion, we have observed 
only the settling of the QSS in its inital form, before the 
slow collisional evolution has caused any significant vari- 
ation in the direction of the BG equilibrium. 




500 



Fig. 22. Magnetization as a function of time M(t) for E — 
0.58. It shows damped oscillations. They may be related to 
Landau damping for inhomogeneous QSS. 



4.1 Gaussian initial conditions 

In this case, for a given energy E, the initial velocities are 
extracted from the distribution 



1 



^/tt (AE - 2; 



: exp 



V 

' AE-2 



(92) 



The critical energy for Vlasov stability is Ec = 0.75. We 
have performed simulations at several energies smaller 
than and precisely at the energy values 0.51, 0.55, 
0.58, 0.60, 0.62, 0.65 and 0.69. We note that the smallest 
possible energy for an initial condition with M = is 0.5 
(see Eq. ([Ts])). In all cases, after the rapid relaxation, the 
system reaches a QSS with a magnetization smaller than 
the equilibrium value, and thus also with an average ki- 
netic temperature smaller than at equilibrium (again from 
Eq. (18 1 one clearly sees that, for a given energy E, a de- 

we 



(91) crease of M implies a decrease of T^m)- In Figure 22 



show a representative plot of the magnetization vs time, 
M{t), for the run at E = 0.58. 

It appears that, after the rapid increase of M and af- 
ter a short transient with large oscillations, the system 
settles to a state with small persistent magnetization oscil- 
lations. The almost regularity of the oscillations indicates 
that they are not due to finite size noise, but that they 
are an intrinsic property of the QSS. Therefore, strictly 
speaking, we should consider the following two possibil- 
ities for the state reached by the system after the rapid 
relaxation. Either the state is not a stationary state of the 
Vlasov equation, but rather a quasi-periodic state; or the 
oscillations are due to a sort of Landau damping of an 
inhomogeneous QSS, and they will eventually die out on 
a time scale large with respect to the one simulated here. 
The first case is reminiscent of the results found by Morita 
and Kaneko [H], that studied initial conditions different 
from ours, and in some cases they obtained almost peri- 
odic oscillations of the magnetization. 

A similar picture holds for the other energies studied, 
in some case with oscillations of a somewhat larger am- 
plitude. Whatever the explanation of the oscillations, we 
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Fig. 23. The distribution /(e) has been fitted by a pure poly- 
trope n = 3/2 (long-dash) or n = 2 (short-dash). We have also 
plotted the initial condition (dot). 
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Fig. 24. The distribution /(e) has been fitted by a polytropic 
core n = 3/2 (long-dash) or n = 2 (short-dash) and a halo. We 
have also plotted the initial condition (dot). 



can use the average value of the magnetization (and con- 
sequently the average value of the kinetic temperature) to 
characterize the state of the system; since the oscillations 
are small we can approximately consider the state to be 
stationary. 

As already noted in the Introduction, a steady state 
(stable or unstable) of the Vlasov equation is characterized 
by a one particle distribution function that depends on 
and V only through the individual energy e = v'^ /2 + <P{0). 
Then, in our QSS the coordinates (6*, v) of the particles 
should be distributed according to this property. In our 
simulations, where My = 0, the potential energy 'P{9) is 



given by 1 — cos (see Eq. ( 12 )), where is the aver- 
age value of the magnetization mentioned above. We have 
therefore divided the one particle phase space in cells small 
enough to characterize each cell with a good approxima- 
tion with a single value of the individual energy e (e.g., 
the value of the central point of the cell) , but big enough 
to contain a large number of particles. Taking a snapshot 
of the system configuration in the QSS (at the end of the 
run), we have then counted the number of particles in 
each cell, and plotted this number vs the energy value at- 
tributed to the cells. If the points of this plot are arranged 
on a continuous line, this means that the one particle dis- 
tribution function depends only on e, confirming that the 
state is a QSS. 



In Figures 23 and 24 we plot the distribution /(e) for 
the two cases E = 0.58 and E = 0.55, respectively. It 
is clear that, with a very good approximation, the points 
are arranged along a single line in both cases, confirming 
that the state of the system is a QSS. In particular, we 
note that the line is practically a straight line in Figure 
[23} while it is a straight line in Figure [24] up to a certain 
energy e, above which it deviates. 

Looking at the polytropic distribution, Eq. (p6t, and 



at the relation between the indices q and n in Eq. (27), 
we see that an /(e) given by a straight line between the 
minimum and maximum values of the individual energy, 
and zero for higher energies, corresponds to a polytrope 
with n — 3/2. In Figure 24 and in Figure [23] the fit of the 
numerical distributions /(e) with a straight line, i.e., with 



a polytrope with n — 3/2, is reported. However, we report 
also the fit with a polytrope with index n — 2, that, in 
spite of the different functional form, is numerically very 
close to the n = 3/2 fit (it is actually slightly better). 
The use of the n — 2 polytropes is in addition justified by 
the numerical caloric curve that we are going to discuss 
shortly. 

From the operative point of view the fit is performed 
by writing Eq. ( 26 1 for g > 1 in the following form: 



/(e)=A(6-e)'; 



-1/2 



(93) 



with n > 1/2. Once chosen, from the numerical distri- 
bution data, the minimum and maximum energies within 
which the fit has to be performed {emin and emax respec- 
tively), the parameter b is equal to emax, while A{b — 

Emin)" IS the valuc of the fit at the minimum energy. 

The difference between E = 0.55 and E = 0.58 is evi- 
dent: for E — 0.55 the part of /(e) at the higher energies 
e is outside the fit. It is due to a halo of particles, and 
we have checked that this halo is robust and does not dis- 
appear at later times. Later, showing the results for the 
waterbag initial conditions, we will present a picture of the 
location of the particles in the one particle phase space in 
another case in which the halo in the distribution function 
/(e) is evident. 

The numerical distributions /(e) for the other energies 
are not plotted here, and we briefly describe what has 
been obtained. For the smallest energy, 0.51, and for the 
energies 0.60 and 0.62, the fit with n = 3/2 or n = 2 is 
equally good, with the exception of the higher energies e 
for E = 0.51, where a halo is present as for E = 0.55. For 
E = 0.62 a small bending at the smaller and at the higher 
energies e begins to appear in the numerical distribution, 
such that its second derivative is negative at small energies 
and positive at high enegies. This bending is marked at 
E — 0.65 and even more pronounced at E = 0.69. For 
these cases, the fit with a polytrope is not good, especially 
for E = 0.69. This can easily be understood from Eq. (93): 
the sign of the second derivative of the polytrope is the 
same in all the energy range [emin, ^max], i-e., positive for 
n > 3/2 and negative for n < 3/2. 
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Fig. 25. The numerical caloric curve displays a region of neg- 
ative kinetic specific heat Ckin = dE/dT^in < 0. It has been 
fitted by a polytrope n — 3/2 for which Ckin = —43/34 and 
by a polytrope n = 2 for which Ckin = —5/2. 



In Figure [25] we plot, with white circles, the kinetic 
temperature as a function of the energy for the QSS 
reached by the system with gaussian initial conditions. 
Accordingly with the fit performed in Figures [23] and [24] 
we fit this numerical caloric curve with the analogous one 
that holds for n — 3/2 and n = 2 polytropes (the fit by 
n = 2 polytropes appears to be very good). Both of them 
show a negative kinetic specific heat in that energy range, 
and so does the numerical caloric curve. The explanation 
of this negative specific heat region in terms of polytropic, 
or close to polytropic distributions, is an important result 
of our paper. This is intrinsically due to incomplete relax- 
ation (lack of ergodicity). We "guess" that the Lynden- 
Bell theory, which assumes complete mixing, would not 
produce a negative specific heat regiorf^ 

We see that at the higher energies, 0.65 and 0.69, the 
fit of the caloric curve is less good, coherently with what 
has been found for the fit of the distribution functions. 
We note that at these energies the kinetic temperature 
of the QSS is not much higher than that of the homo- 
geneous branch (equivalently, the magnetization of the 
QSS is small). It is then plausible that the relaxation from 
the homogeneous initial condition to the stable state can 
be approximately described by a quasi-linear theory, that 
treats the inhomogeneity of the distribution function as a 
perturbation of the angle-averaged homogeneous distribu- 
tion. This is left for a future investigation [59j . 



4.2 Semi-ellipse initial conditions 

Now, for a given energy E, the initial velocities are ex- 
tracted from the distribution 



1 



n {4:E - 2) 



\/8E -4-1)2 



(94) 



^® We have not computed the prediction of the Lynden-Bell 
theory which would require to develop a specific algorithm to 
treat multi-levels initial distributions. However, the Lynden- 
Bell distribution has not a compact support so it does not 
account for observations. 
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Fig. 26. The distribution /(e) has been fitted by a pure poly- 
trope n — 1 (dash). We have also plotted the initial condition 
(dot). 
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Fig. 27. The distribution /(e) has been fitted by a polytropic 
core n = 1 (dash) and a halo. We have also plotted the initial 
condition (dot). 



corresponding to a n = 1 homogeneous polytrope. The 
critical energy for Vlasov stability is i5c = 5/8 = 0.625. 
We have performed runs at the energies 0.51, 0.55, 0.58, 

0. 59, 0.60 and = 0.61. As before, after a rapid relaxation 
the system reaches a QSS with a magnetization and a ki- 
netic temperature smaller than the respective equilibrium 
values. 

We do not show analogous plots of the magnetization 
vs time, since they are very similar to those obtained for 
the gaussian initial conditions, with small persistent os- 
cillations. We go directly to the plots of the numerical 
distribution functions. We show the two cases E — 0.59 
and E — 0.55, respectively in Figures [26| and [27] 

The points of the numerical distribution functions are 
arranged along a line, as it should be for a QSS. As for the 
gaussian case, the distributions at the smaller energies E 
present a halo. The fit with a polytropic distribution has 
been performed with the n = 1 polytrope. The fit appears 
rather good in both cases, of course except for the halo at 
E = 0.55. The picture is similar as for the gaussian case, 

1. e., there is a halo at small energies (we have one also 
for E — 0.51) and the fit with the polytrope worsen at 
high energies, here E — 0.60 and E = 0.61, when the sign 
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Fig. 28. The numerical caloric curve displays a region of neg- 
ative kinetic specific heat Ckin = dE/dTkin < 0. It has been 
fitted by a polytrope n — 1 for which Ckin ~ ^1/2. 



of the second derivative of the numerical distributions is 
negative at small energies e and positive at high energies. 
Actually, we found that also at energy E — 0.58 the fit 
with a polytrope distribution is not satisfying. 

In Figure [28] we show the numerical kinetic caloric 
curve of the QSS reached by the system, fitted by the 
n = 1 polytropic kinetic caloric curve, that has a negative 
kinetic specific heat in that energy range. We see that the 
polytropic fit is relatively good for most energies. As for 
the gaussian case, at the highest energies the fit is less 
good, and probably the QSS in that case could be better 
explained by a perturbative analysis. 

We remark that the polytropic index that fits the QSS 
is the same as that of the initial velocity distribution, i.e. 
n = 1. 



4.3 Waterbag initial conditions 

For this class of initial conditions, the initial velocities at 
a given E are extracted from the distribution 

/o(.) = ^^^a(V6¥33-H) , (95) 

corresponding to a n = 1/2 polytrope. The critical en- 
ergy for Vlasov stability is Ec = 7/12. We have performed 
runs at the energies 0.51, 0.55, and 0.582. The last ener- 
gies is only very slightly smaller than the critical energy. 
While sharing the previous general picture, i.e., the rapid 
approach to a magnetized QSS, there are some new fea- 
tures. We begin by showng the plots of the magnetization 
vs time for the two energies E = 0.582 and E = 0.55, 
respectively in Figures [29] and [30] 

We note that the magnetization oscillations are more 
marked than previously, becoming very pronounced at 
E = 0.582. This is relatively close to the situation reported 
by Morita & Kaneko [48 . We noted before that one possi- 
ble explanation of the oscillations is a Landau damping as- 
sociated to the inhomogeneous QSS. Now, these large os- 
cillations suggest that this damping could be very weak. It 




Fig. 29. Magnetization as a function of time M (t) for E — 
0.582. It shows undamped (or slowly damped) oscillations. This 
may be related to the absence of Landau damping for the QSS 
(= inhomogeneous waterbag) as for the homogeneous waterbag 
distribution. 
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Fig. 30. Magnetization as a function of time M{t) for E = 
0.55. The amplitude of the oscillations is smaller than in Figure 
|29| This may be due to the presence of the halo. 



is suggestive to interpret this as a continuation with what 
happens for the homogeneous waterbag QSS, for which 
we know that Landau damping is absent (this being due 
to the singularity of the waterbag distribution, that, con- 
trary to distributions strictly decreasing for increasing \v\, 
admits purely real proper frequencies). 

Nevertheless, we proceed as before and plot the numer- 
ical distributions as a function of the individual energy, 
using the average value of M to define the last quantity. 
In Figures [31] and [32] we plot the numerical distributions 
for the cases E — 0.582 and E = 0.55, together with the 
fit with the n = 1/2 polytrope. At £: = 0.582, the QSS is 
an almost pure inhomogeneous waterbag distribution and 
the phase portrait looks like Fig. [21] At the smaller energy 
E — 0.55 there is a halo that cannot be reproduced by the 
polytropic fit. This halo is clearly evident in Figure[33|that 
shows the location of all the particle in the one particle 
phase space in the QSS state a,t E = 0.55. The boundary 
between the more densely populated region and the less 
densely populated region (the halo) is at the energy where 
the halo begins in Figure l32| 
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Fig. 31. The distribution /(e) has been fitted by a "pure" 
homogeneous distribution (waterbag) corresponding to a n = 
1/2 polytrope. We have also plotted the initial homogeneous 
waterbag distribution (dot). 
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Fig. 32. The distribution /(e) has been fitted by a n = 1/2 
polytrope (core) with a halo. This "core-halo" structure is sim- 
ilar to the one found in |36) . We have also plotted the initial 
homogeneous waterbag distribution (dot). 
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Fig. 33. Occupation of the one particle phase space in the 
QSS state at i5 = 0.55. The halo, i.e., the region less densely 
populated, is clearly visible. 
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Fig. 34. Relation between the energy E and the distribution 
value / of a waterbag distribution that is a steady state of 
the Vlasov equation [J^- Homogeneous and inhomogeneous 
distributions are represented. In the simulations, we start from 
an unstable homogeneous waterbag distribution with energy E 
and distribution value fo- We see that only a,t E = Et does the 
inhomogeneous polytrope have the same / as the homogeneous 
polytrope. In the other cases, a halo is necessary to maintain 
the initial value fo in the core of the QSS (as observed in 
the simulations) while conserving energy. However, close to 
E — Et, the effect of the halo is weak and a pure waterbag 
distribution provides a good approximation of the QSS. 



As for the semi-ellipse initial conditions, we remark 
that the polytropic index of the QSS is the same as that 
of the initial velocity distribution, i.e. here n = 1/2. A 
difference with respect to the two previous cases is that 
the fit appears to be good also very close to the critical 
energy Ec- Furthermore, the value of the numerical distri- 
bution function in the QSS is practically the same as that 
of the initial waterbag distribution (see Figures 31 and 
32 1. This means that the core does not mix at all. 



Using 

tms observation, we can explain the presence or the ab- 



sence of the halo. To that purpose, we plot in Fig. 34 the 
relation between the uniform distribution value / and the 
energy E of a pure waterbag distribution that is solution 
of the Vlasov equation (the construction of this Figure is 
explained in [37] ) . This Figure shows the first order phase 
transition between homogeneous and inhomogeneous wa- 
terbag distributions discussed in Sec. |3.8| We see that, in 



general, the homogeneous and inhomogeneous waterbag 
distributions with the same f have a different energy (see 
the vertical line in Fig. 34). Therefore, if / is the same 



in the initial homogeneous waterbag distribution and in 
the inhomogeneous waterbag QSS (as it turns out to be), 
there must necessarily exist a halo of particles in order 
to satisfy the conservation of energy. The halo should be 
particularly important at low energies E (as in Fig. 32 for 
E = 0.55) where AE is large. By contrast, at the transi- 
tion energy Et, the homogeneous and inhomogeneous wa- 
terbag distributions have the same / and E so that the 
presence of a halo is not required. By continuity, close to 
the transition energy (as in Figure 31 for E ~ 0.582), the 
halo should be modest since AE is small (for E — 0.582 
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Fig. 35. The numerical caloric curve has been fitted by a poly- 
trope n — 1/2. Strictly speaking, this caloric curve exhibits a 
first order phase transition at Et = 0.588, leading to Ckin = 0. 
However, metastable states are equally relevant as fully stable 
states, so the whole series of equilibria has been plotted. 



we find AE/E ~ 0.01). These arguments are consistent 
with the observations. 

In Figure [35] we plot the numerical kinetic caloric 
curve, together with the n — 1/2 kinetic caloric curve. In 
principle, the n = 1/2 kinetic caloric curve should display 
a first order phase transition as shown in Fig.[T5j However, 
for systems with long-range interactions, the metastable 
states (local entropy maxima) are extremely robust, and 
they are as much relevant as fully stable states (global 
entropy maxima) . For that reason, we have chosen to rep- 
resent the full series of equilibria, displaying both global 
entropy maxima, local entropy maxima, and even saddle 
points of entropy. 

The first thing to note is that, apart from the highest 
energies, the kinetic temperature is very close to that of 
the BG equilibrium. However, as clearly proved from the 
numerical distribution functions, the state is far from be- 
ing the BG one, that has a Boltzmann distribution. The 
second thing to note it that, also in this case, there is a 
region of negative kinetic specific heat. 



A striking feature of Figure 32 is the core-halo state. 
This core-halo state, arising from a waterbag initial con- 
dition, was previously observed by Pakter & Levin |36l , 
although they did not explicitly calculate the curve /(e) 
(they observed the core-halo state from the phase space 
portrait). A new contribution of our work is to show that 
this core-halo state is also present at l ow e nerg ies fo r other 
types of initial conditions (see Sees. 4.1 and 4.2). In all 



the cases considered, the core can be fitted by a polytrope 
with a different index (n — 1/2 for the waterbag initial 
condition). This generalizes the results of [36] . 

There are a few differences between our approach and 
the approach of [SS]. First, we have considered a waterbag 
initial condition with Mq = while they took Mq = 0.40. 
For Mq = 0.40, the Lynden-Bell theory predicts a second 
order phase transition (see [19] or Figures 2 and 4 of 211) 
which is in clear disagreement with the first order phase 
transition reported in [35]. By contrast, for Mq — 0, the 
Lynden-Bell theory predicts a first order phase transition. 
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Fig. 36. Full line: The distribution /(e) obtained in the QSS 
reached by the waterbag initial condition with E — 0.6 and 
Mo = 0.1711 (we have indicated by a dotted line the value 
of the initial distribution function /o). Dashed line: the distri- 
bution predicted by the Lynden-Bell theory. We see that the 
agreement is almost perfect. We may also plot ln[//(/o — /)] 
as a function of e, as suggested in [2^, and check that it is a 
straight line. In this way, we do not have to compute the the- 
oretical Lynden-Bell distribution. The Lagrange multipliers a 
and P may be determined from the relation ln[//(/o — /)] — 



It is, however, different from the first order phase transi- 
tion reported in Figure [35] because the Lynden-Bell predic- 
tion corresponds to a distribution that is partially degen- 
erate while we find a distribution that is either completely 
degenerate, or with a core-halo structure. Therefore, both 
in |36| and in the present study, the Lynden-Bell predic- 
tion fails although the situation is a bit different. Secondly, 
in order to obtain their theoretical magnetization curve, 
Pakter & Levin f36J assume that the distribution function 
in the core of the QSS is equal to the initial distribution /g 
(our numerical distribution functions give further support 
to this assumption, as shown in Figures 31 and 32) and 
determine the properties of the halo by a semi-analytical 
approach. While we agree with their procedure which pro- 
vides a good prediction of the QSS for all energies, we have 
proceeded differently. We have obtained the theoretical 
magnetization curve (or kinetic caloric curve) by assum- 
ing that the QSS is a pure polytrope n = 1/2. Since we 
ignore the halo, the distribution function fth that we theo- 
retically compute is generally different from the initial dis- 
tribution function /q in order to satisfy the conservation of 
energy (see the horizontal line in Fig. 34 ) . This procedure 
provides a reasonable agreement with the numerical sim- 
ulations in the region where the halo is not pronounced, 
i.e. close to the transition energy, where Af / fg ^ 1 (for 
E = 0.582 we find Af/f = 0.02). By contrast, it clearly 
fails for lower energies where the halo is significant. This 
simply reflects the fact that a pure polytrope cannot hold 
for all the energies, and that a core-halo state is required. 

In this paper, we have chosen to focus on examples 
where the Lynden-Bell prediction fails. However, as re- 
called in the Introduction, we know that the Lynden-Bell 
predictions are in many cases verified. In order to mod- 
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erate our message about the inadequacy of the Lynden- 
Bell distribution to describe the QSS in certain cases, 
we give here an example of a QSS in which the numer- 
ical distribution function agrees with that predicted by 
the Lynden-Bell theory with extremely good precision. In 
Figure [36] we show the numerical distribution function in 
the QSS reached by the system initially prepared with a 
waterbag distribution of the same kind as those consid- 
ered in [23|; precisely, we considered a rectangular wa- 
terbag initial distribution at E = 0.6 and initial magne- 
tization Mq — 0.1711 (this initial condition corresponds 
to the point E — 0.6 and /o — 0.113 on Figure 10 of 
[3T]). The Lynden-Bell theory predicts in this case a mag- 
netized QSS. Not only do we find the magnetization in 
the simulation agrees with that of the Lynden-Bell the- 
ory {Mqss = 0.348), as observed previously in [IJj, but 
we also find that the numerical /(e) is that predicted by 
the theory (a feature that was not explicitly checked in 
[H]). The Lynden-Bell function is plotted, in Figure 36 
together with the numerical distribution function. 



5 The approach to BG equilibium of 
homogeneous Vlasov stable states 



In previous works [35)57j . we have studied the modalities 
of the approach to BG equilibrium of the system prepared 
in a Vlasov stable homogeneous state at energies below the 
thermodynamical critical energy Ec = 0.75. In that case, 
we were interested in the lifetime of these homogeneous 
states, whose slow evolution is governed by the finite size 
effects. We know that this evolution changes slowly the 
state of the system, that remains homogeneous, until the 
distribution becomes Vlasov unstable and the system be- 
gins to magnetize and to approach equilibrium. In Ref. 
[57] we showed that, during the slow "collisional" evo- 
lution, and for homogeneous distributions, the velocity 
distribution f{v,t) can be fitted with good approxima- 
tion by a polytropic function, whose index n(t) changes 
with time in correspondence with the change of the dis- 
tribution. We considered in particular the mostly studied 
energy E = 0.69, preparing the system with velocities ex- 
tracted from a semi-elliptical n = 1 polytropj^ This ini- 
tial state is Vlasov stable, since E > Ec — 5/8. We found 
that, during the dynamics, the index n{t) of the fitting 
polytrope increases. Solving Eq. ([90|) for n, we see that 



for a given energy E < 3/A the homogeneous polytrope is 



^"^ In [35], the authors started from a rectangular waterbag 
distribution with energy E — 0.69 and vanishing magneti- 
zation Mo = 0. This initial condition is Vlasov stable since 
E > Ec = 7/12 so it does not experience phase mixing and 
violent relaxation. However, it slowly evolves due to finite N 
efi^ects. They found that the system rapidly forms a velocity 
distribution with a semi-elliptical shape. This can be inter- 
preted as a polytrope n = 1 [29]. Yamaguchi et al. [54] had 
previously obtained the same result but they did not to recog- 
nize the polytrope (Tsallis distribution) ; see discussion in [22 . 
In Ref. [37], we directly started from the polytrope n = 1 to 
accelerate the simulation. 




3x105 



Fig. 37. Magnetization vs time for homogeneous initial con- 
ditions with semi-elliptical velocity distribution, at -E = 0.69, 
for a system with A*' — 2^"^ particles. 



Vlasov stable if, and only, if n < nc = (4i? - 2)/(3 - AE). 
For E = 0.69 this gives n < — 19/6. We found that 
when n approaches the critical value ric the homogeneous 
distribution begins to destabilize, and the system begins 
to magnetize. 

In Ref. |57j we had not analyzed the distribution func- 
tions after the homogeneous phase becomes Vlasov unsta- 
ble, and the system becomes magnetized. Here, we are in- 
terested to know whether the magnetized states /(e, t) can 
still be fitted by inhomogeneous polytropes with a time 
dependent index n{t). If true, the index will start from 
about ric = 19/6 and increases towards infinity, which 
corresponds to the BG state, as time goes on. From the 
study of (lOj , we know that the inhomogeneous polytropes 
with index n > umcp — 0.68 are always Vlasov stable 
(see Remark 1). Therefore, we conclude that the whole 
sequence of inhomogeneous polytropes with index larger 
than Tic = 19/6 is Vlasov stable. 

Here, we have performed a simulation of the system ini- 
tially prepared, at E = 0.69, in a homogeneous state with 
the velocities distributed according to the semi-elliptical 
distribution. This is the same initial condition considered 
in |57| . but now we have studied the one-particle distri- 
bution function /(e, t) in the time range mentioned before 
(i.e., in the magnetized phase). We could not use the same 
number of particles as for the analysis of the QSS reached 
from a Vlasov unstable state (see Sec. [4]) since this would 
have required an unmanageable computer time. Our runs 
have been made with 2^^ particles, representative of a sys- 
tem of 2^^ particles, for the same symmetry property as 
before. 

In Figure[37j we plot the time evolution of the magneti- 
zation. We remark that at the end of the run the BG state 
has not yet been reached (the equilibrium magnetization 
is about 0.31), but the length is sufficient for our analysis. 
Here, we consider the one particle distribution functions 
during the time range in which the magnetization rises. 

We have computed the distributions at the times t = 
1.5 X 10^ t = 1.7 X 10^^= 2.0 X 10^ t = 2.4 X 10^ t = 
2.8 X 10^. From Figure 37 it can be seen that these times 
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Fig. 38. One particle distribution function at time t = 2.8 x pjg^ Magnetization vs time for waterbag initial conditions 
10 . It has been fitted with a polytrope with index n = 4.3. with A/(0) = 1 at i5 = 61 



38 



we 



are all after the destabilization of the homogeneous phase 
(the corresponding values of the magnetization, averaged 
over a short interval of time, are M = 0.17, M — 0.19, 
M = 0.26, M = 0.26, and M = 0.28). In Figure 
show the distribution /(e) obtained sX t = 2.8 x 10 

In all cases, the points are arranged along a line and 
therefore at all times the system is in a stationary state 
of the Vlasov equation. Therefore, for all the time range 
during the approach to BG equilibrium, the system passes 
by a sequence of quasi stationary steady solutions of the 
Vlasov equation, slowly evolving under the effect of "col- 
lisions". This property is due to the scale separation be- 
tween the relaxation time (larger than N) and the dynam- 
ical time (of order unity). Furthermore, the system never 
goes through Vlasov instability, and all the dynamics is 
governed by the coUisional finite size effects. 

The distributions have been fitted with polytropic 
functions, with index n = 2.5, 3.0, 3.5, 4.0 and 4.3, re- 
spectively. The fit is very good in all cases (it is shown in 
Figure 38 for the case n = 4.3). In addition, the fits sug- 



gest an increasing index with time. This is of course nat- 
ural since the BG equilibrium corresponds to a polytrope 
with index n = oo. We note that the first index calculated 
in the magnetized phase is somewhat smaller than ric , and 
the second index is very close but slightly smaller than ric. 
In principle this can be considered perfectly legitimate; in 
fact, even though the numerical distribution function is all 
the time very close to a polytrope, small discrepancies can 
give rise to small uncertainties in the value of the index 
n, considering that in the case of inhomogeneous distribu- 
tion we have to fit a function /(e) instead of a function 
f{v), and there are more sources of numerical errors (like, 
e.g., the determination of the magnetization M in the def- 
inition of the individual energy e). Furthermore, a slight 
change of the maximum energy of the fit could correspond 
to a slight change of the polytropic index n. We remark 
that for small values of n, like those of the previous Sec- 
tion, this problem of the small uncertainty in the value of 
n is much less relevant, because of the steeper decrease of 
those polytropes. It is also possible that a value of n{t) 
smaller than ric in the inhomogeneous phase is real (i.e. it 
is not a artifact due to the fit). We may well imagine that 



the homogeneous phase destabilizes when n{t) — ric and 
that the inhomogeneous polytropes just after the transi- 
tion have an index smaller than ric. In other words, the 
evolution of n{t) could be non-monotonic. 



6 Waterbag initial condition with Mq = 1 

In this Section, we consider a different class of initial con- 
ditions, i.e., a waterbag distribution for the velocities, and 
all the angles at 6* = so that the initial magnetization 
is Mq = 1 (this distribution is unsteady). That was the 
first initial condition considered for the HMF model, the 
one that revealed the existence of QSSs [8jll| 12j. These 
conditions have been extensively studied, with the pur- 
pose to determine the scaling with N of the QSS lifetime 
and of its magnetization. It is known that there are large 
fiuctuations from run to run, so it is necessary to per- 
form an average on several runs to estimate the mentioned 
scaling. In Ref. it was shown that the fluctuations 
can be reduced by using the so-called isotropic waterbag 
conditions, in which the velocities are not randomly ex- 
tracted from the waterbag, but are taken equally spaced. 
This is equivalent to using normal waterbag distributions 
and performing averages over many runs. We adopt this 
strategy here, but now with the purpose to study the one 
particle distribution functions in the QSS. 

We have performed runs at the energies E = 0.55, 
E ^ 0.57, E = 0.61, E = 0.64, E = 0.65, E = 0.66, 
E = 0.665, E = 0.67 and E = 0.69. As in Section [4) we 
simulate a system with 2^® particles. The last energy is 
the one mostly studied in the literature, where it is shown 
that the magnetization of the QSS vanishes in the large N 
limit. We confirm this result here. On the other hand, we 
found that at all the smaller energies the magnetization of 



the QSS does not vanish. In Figures 39 and 40 we show the 
time evolution of the magnetization for the two energies 
E = 0.61 and E — 0.69, respectively. 

In the early stages of the QSS, the occupation of the 
one particle phase space presents holes that tend to dis- 
appear. This peculiarity is shown in some figures where 
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Fig. 40. Magnetization vs time for waterbag initial conditions 
with M(0) = 1, at _B = 0.69. 
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Fig. 41. Occupation of the one particle phase space in the 
early stage {t = 100) of the QSS at 
waterbag initial conditions for the vi 
The plot shows a hole in the center 
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Fig. 42. Occupation of the one particle phase space at later 
times (t = 500) of the QSS at E = 0.61. The hole at the center 
has shrunk considerably. 
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early stage {t = 100) of the QSS at E = 0.61, for the run with Fig. 43. Distribution function /(e) in the QSS for E = 0.61. 
waterbag initial conditions for the velocities, and with Mo = 1. It has been fitted, neglecting the high peak at small energies, 

by a no — 3.56 polytrope. 



we plot the location of all the particles in the one parti- 
cle phase space. The plots correspond to the two energies 
E = 0.61 and E = 0.69, and to two different times for 
each energy. Figure [41] represents the phase space occu- 
pancy a.t E = 0.61 at the time t = 100, which is one fifth 
of the duration time of the run, as it can be seen from 



Figure [39j The phase space occupancy at the latest time, 
i = 500, is plotted in Figure |42| 

Although the hole at the center of the phase space 
has shrunk considerably a,t t = 500, it is still very clearly 
visible. This appears to be in apparent contrast with the 
numerical distribution function /(e) at the same time t = 
500, that is shown in Figure |43j As in the previous cases, 
the disposition along a line of the numerical data confirms 
that the state is a QSS. However, the interesting point to 
stress is the presence of a pronounced peak at the smallest 
energies e. From the expression of the individual energy 
e, we know that the smallest energies correspond to the 
particles located near the center of the one particle phase 



space. Then, we could argue that there is a contradiction 
between this high peak in /(e) and the hole in the center 
of the phase space. 

To clarify this point we have enlarged the region near 
the center of the phase space (we do not report here these 
plots), and this has revealed that the hole contains in fact 
many particles arranged along thick filaments, that in the 
two-dimensional plots embracing the whole phase space 
are not visible. These filaments could be explained by the 
following argument. 

We know that, as long as the evolution of the one 
particle distribution function f(9,v,t) is governed by the 
Vlasov equation, the mass of each of its phase levels is 
conserved. However, the evolution causes a complex inter- 
twinement of the various phase levels; this is what allows 
to consider coarse-grained distributions, and justifies the 
fact that the numerical distributions, that inevitably have 
an implicit averaging (at small scales as it may be) have 
to be compared with the theoretical coarse-grained distri- 
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Fig. 44. Occupation of the one particle phase space in the 
early stage {t = 100) of the QSS at _E = 0.69, for the run with 
waterbag initial conditions for the velocities, and with Mq = 1. 
The plot shows two holes. 



V 



t t 



Fig. 45. Occupation of the one particle phase space at later 
times {t = 500) of the QSS at £ = 0.69. The holes have disap- 
peared. 
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Fig. 46. Distribution function in the QSS for E = 0.69. It 

has been fitted, neglecting the high peak at small energy, by a 
n — 1 polytrope. 




Fig. 47. Numerical caloric curve for the QSS reached from the 
waterbag AI{0) = 1 initial conditions. It has been fitted by a 
polytrope with index no — 3.56 for which Ckin — >■ oo close to 
the bifurcation point. The fit is relatively good except for three 
points that reveal a sort of "hump". This hump seems to be 
real (i.e. it is not a numerical artifact). It presents a region of 
negative kinetic specific heat. 



buttons. For the case under study in this Section, the fine- 
grained distribution function has only two levels, zero and 
infinite (i.e., it is a (5 function), since all the particles are 
initially at = 0. It may then be argued that for a while 
there remains very high values also for the coarse-grained 
distribution function, that are gradually decreased by the 
mixing. Unless the coarse-grained distribution is analyzed 
at a very small scale (but this would be possible only with 
a very large number of particles) , high values can be ob- 
tained only with particles arranged approximately along 
a one-dimensional region, i.e., a line. 



In Figures 44 and 45 we show the occupancy of the 
one particle phase space for the energy E — 0.69, at the 
times t = 100 and t — 500 respectively. Again, at the 
earlier time the plot presents holes, although they are not 
at the center. We remind, however, that now the QSS is 



not magnetized, as shown in Figure 40 At the later time 
the holes have disappeared. 



Also in this case the one particle distribution function 
a.t t — 500 has a pronounced peak at the smaller energy 
e, as it can be seen in Figure 46 As a matter of fact, the 
peak is present for all energies E that have been simulated. 
However, neglecting the high peak at small energy, we 
see that the distribution function at i? = 0.69 can be 
fitted with a polytrope with n = 1. In a previous work 
[35) , with numerical simulations performed with much less 
particles, and therefore with stronger finite size effects and 
faster QSS evolution, the peak probably had disappeared 
very soon, and no evidence of it was found (the velocity 
distribution f{v) could be nicely fitted by a pure n — I 
polytrope) . 

It is interesting to plot the numerical kinetic caloric 
curve obtained from our simulations. This is done in Fig- 
ure 47 Except for the curious points around the energy 
E — 0.66 (more precisely, the points at i5 = 0.65, E = 0.66 
and E ~ 0.665) that reveal a sort of "hump", it seems 
that Tkin is almost constant with E. This "plateau" cor- 
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responds to a region with infinite kinetic specific heat 
Ckin = oo. Actuahy, the numerical caloric curve is well- 
fitted by the physical caloric curve of no ~ 3.56 polytropes 
that presents a region of infin ite s pecific heat close to the 
bifurcation point (see Section 3.6). The distribution func- 



tion of the QSS at E = 0.61 is also well-fitted by a uq 
polytrope, except again for small energies (see Figure 43). 
These new results do not support the "guess" made in our 
previous paper |10j that the numerical caloric cuve can be 
fitted by n = 1 polytropes. They arc also in conflict with 
the early claims of |8lllj that the numerical caloric curve 
presents a region of negative kinetic specific heat (this re- 
gion is not clearly apparent in our simulations, except in 
the hump). However, we know that a system prepared in 
the waterbag distribution has a strange dynamics, so that 
more averages may be needed to obtain reliable results. 



7 Conclusion 

In spite of the large amount of numerical studies per- 
formed on the HMF model, there are two main prob- 
lems that are still poorly understood. The first one is the 
anomalous scaling law with N of the homogeneous QS^^ 
the second, treated in this work, is the dependence of the 
QSS reached by the system from a generic initial condi- 
tion. The first problem should be restricted to the ho- 
mogeneous QSSs of one-dimensional systems, where the 
coUisional effects vanish at order The second prob- 

lem is expected to be relevant to long-range systems in 
general. 

The Lynden-Bell theory determines, at least in prin- 
ciple, the distribution function of the QSS, although a 
technical difficulty is to compute it in the generic case of 
multi-levels initial distributions. On the other hand, var- 
ious numerical simulations have shown that the theory 
works only in some cases, failing in others. Since the the- 
ory is based on ergodicity or efficient mixing, the reason of 
the failure has to be traced to an incomplete mixing |28] . 
If the system does not mix well, we should not expect any 
"universality": The QSS will depend on the degree of mix- 
ing. This rises the level of difficulty of the problem since 
we do not have indications on the further constraints that 
are imposed by the lack of full mixing. 

In this work we have not offered a solution of this 
problem, in the sense that we have not proposed a way 
to predict the QSS as a function of the initial condition. 
Rather, we have proposed to analyze the correspondence 
between a class of stationary states of the Vlasov equation 
and the numerical QSS reached by the system starting 
from various kinds of initial conditions. We have found 
numerically that, in many cases, the QSS is well-fitted 
by a polytropic distribution with a compact support. The 
fact that the distribution has a compact support is con- 
sistent with the phenomenology of incomplete relaxation 
according to which some parts of the phase space (usually 
those with high energy e) are not sampled by the system. 



In this respect, the value of the index n may be a measure 
of mixing. 

Polytropic distributions were introduced in different 
situations: 

(i) In astrophysics, they were introduced as particular 
steady states of the Vlasov equation [33]. Furthermore, it 
was shown that they extremize a "pseudo-entropy" func- 
tional S = - J /i+i/("-3/2) ^j.^^ at fixed mass M and 
energy E |60| , and that maxima of such functionals are 
dynamically Vlasov stable [61l62j . This is the interpreta- 
tion that we have followed in this paper. 

(ii) Polytropic distributions were also introduced by 
Tsallis '13| from an approach of generalized thermody- 
namics. In that case, the functional S is interpreted as a 
generalized entropy and its maximization at fixed mass 
and energy determines the most probable state of a non- 
ergodic system for which some regions of phase space are 
forbidden by the dynamics. Several classes of generalized 
entropies are possible but the Tsallis entropy is the sim- 
plest deviation to the Boltzmann entropy. We leave open 
the debate about the relevance of this interpretation and 
refer to |31I52] for more discussions about this issue. 

In this work, we have performed numerical simulations 
of the HMF model initially prepared in a Vlasov unsta- 
ble steady state, and we have analyzed the properties of 
the QSS reached by the system after the rapid relaxation. 
There have been previous studies dedicated to the descrip- 
tion of the QSS in the HMF model, and the following ones 
have treated the subject under different conditions. In 
Refs. |19I21I23I26| . the authors have compared the numeri- 
cally obtained QSS phase diagram with the theoretical one 
arising from the application of the Lynden-Bell theory to 
rectangular waterbag initial conditions. They found that 
the results agree in some regions of the phase diagram, 
and disagree in others. Morita and Kaneko [IH] have ana- 
lyzed the magnetization of the QSS reached by particular 
initial conditions, i.e. of the same form as that attained in 
the BG equilibrium, but with out-of-equilibrium param- 
eter values; the value and the possible oscillations of the 
QSS magnetization have been studied as a function of the 
parameters of the initial distribution. Finally, Pakter and 
Levin |36| have studied the same two levels initial distribu- 
tions as in |19| , with the purpose to propose an alternative 
QSS one particle distribution function, with respect to the 
Lynden-Bell functions, for the cases where the latter does 
not work. 

From the point of view of the initial conditions ana- 
lyzed, our work is in a sense complementary to those just 
cited. We have studied several classes of initial distribu- 
tion functions, some of which with only two levels, as in 
19 21 23 26 and [36 , and some with continuous values, 
as in |48) . In addition we have analyzed the properties 
of the QSS distribution functions, showing that they can 
be fitted in many cases with polytropic functionj^ In 
this respect, the spirit of our work is closer to [36" , where 
the authors try to find the expression for the distribution 



Some interesting results in this direction have been ob- 
tained recently [68] . 



Similarly, some QSSs in 2D turbulence have also been fitted 
by polytropic (Tsallis) distributions [63164165] , although this fit 
cannot be universal 1661. 
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function. Our results for the QSS reached from waterbag 
homogeneous initial conditions confirm the results of Ref. 
|36| : we found that the QSS distribution function has a 
core-halo structure; the limits of the core, in the one par- 
ticle phase space, are given by a constant energy line, and 
the core can be fitted by a function /(e) with a waterbag 
structure, i.e., a n = 1/2 polytrope. In addition, we found 
that the magnetization in the QSS presents oscillations, 
reminiscent of those found in [48] . This could be due to the 
absence of Landau damping, analogously to what holds 
for the homogeneous waterbag function, characterized by 
purely real proper frequencies for the corresponding lin- 
earized Vlasov equation. 

We found that polytropic functions can fit the QSS dis- 
tributions also in the other cases, with an index depending 
on the class of initial condition, although the fit is not very 
good for all the energies considered. In particular, for en- 
ergies slightly smaller than the instability threshold the 
fit appears less satisfying, and probably those cases could 
be considered with a quasilinear approach [55] ■ The fit 
with polytropic functions has allowed us to compare the 
numerical caloric curves with the kinetic caloric curves 
computed on the basis on the theory of polytropes [lOJ. 
In particular, using this theory, we have been able to ex- 
plain the negative kinetic specific heat region observed in 
the numerical simulations for intermediate energies. This 
is one of the main results of the paper. For lower energies, 
the system is not a pure polytrope but it takes a core-halo 
structure. The core can be fitted by a polytrope, which 
may differ from the waterbag distribution. The halo may 
also be fitted by a polytrope so that the core-halo state 
may be viewed as a "mixture" of two polytropej^ There- 
fore, our results complement, and generalize, the findings 
in Ref. 36J for the waterbag initial condition. 

An interesting problem is to understand why the QSS 
can be fitted by a core-halo state with a polytropic core. 
Furthermore, the polytropic index of the QSS seems to 
be correlated to the initial distribution function. Indeed, 
for the semi-ellipse and waterbag initial conditions, we 
find that the QSS can be fitted by a polytrope with index 
n — 1 and n — 1/2, respectively. This could be interpreted 
by saying that the system keeps memory of the initial con- 
dition. For the waterbag case, the core does not mix well 
so it keeps its uniform distributiorj^ On the other hand, 
a halo is formed probably due to parametric resonance 



For example, in Fig. 
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the distribution function in the 
halo is approximately constant so that the core-halo state may 
be fitted by two polytropes n = 1/2 with different valu es o f / 



We remark that, in all considered cases (see Figs. 
261 [271 [STI and 
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32 1, the maximum value of the distribution 



function is conserved. This is not a trivial result since we expect 
the coarse-grained distribution function f,^^^ to decrease as 
the system mixes (this is what the Lynden-Bell theory usually 
predicts). This shows that the core does not mix well. This 
property is also observed in self-gravitating systems |46| and 
in two-dimensional turbulence |69| . The conservation of the 
maximum value of the distribution function can be explained 
qualitatively in terms of a kinetic theory (see Appendix A of 
170171) . Appendix B of f^S^, and f30!). 



[36) . As explained in this paper, the halo is necessary in 
order to conserve the value of /q in the core while keep- 
ing the total energy unchanged. Similar explanations may 
hold for more general distributions. 

We have also studied the approach to BG equilibrium 
of a system initially prepared in a homogeneous Vlasov 
stable state. This approach is due to finite size effects. In 
Ref. i57| we found that the slowly changing homogeneous 
distribution function can be fitted by a sequence of poly- 
tropes with a time dependent index. Here we extended 
that work, and we have found that this fit is good also in 
the successive time range, during which the system pro- 
gressively magnetizes. We found that the polytropic index 
increases with time, although the index in the early stages 
of the magnetized phase appears to be slightly less than 
the critical index for the homogeneous polytrope. This in- 
crease is to be expected since in the BG equilibrium state 
the index should approach infinity. The reason why poly- 
tropic distributions fit well the collisional relaxation of 
the system remains a mystery. We recall, however, that a 
similar observation has been made by Taruya & Sakagami 
[57] for self-gravitating systems. Therefore, polytropic dis- 
tributions seem to be ubiquitous in long-range systems. 

Finally, we have considered the QSS reached by a pe- 
culiar initial state, in which all the rotators are at = 0, 
i.e., with the initial magnetization equal to 1. In this case, 
the dynamics is complex due to the formation of phase 
space holes and thick filaments, and the results have been 
more difficult to interpret. We found that, except for a few 
energy values, the numerical caloric curve presents a hori- 
zontal segment corresponding to an infinite kinetic specific 
heat. However, we do not have an explanation for the odd 
results for the three energies that do not comply with the 
others and seem to reveal a hump. We also found that the 
fit with polytropes is here more problematic. For the only 
energy with a non magnetized QSS, E — 0.69, we found 
that the QSS distribution function can be fitted with a 
semi-ellipse, i.e., an n = 1 polytrope, only neglecting the 
pronounced central peak at the smaller energies. We have 
shown that this peak tends to disappear during the slow 
QSS evolution. For the others energies with a magnetized 
QSS, the numerical caloric curve has been compared to the 
kinetic caloric curve of no — 3.56 polytropes that presents 
an infinite kinetic specific heat Ckin = +oo (plateau) close 
to the bifurcation point. The fit is fairly good, except for 
the three anomalous points mentioned above. The distri- 
bution function at E = 0.61 has also been fitted by a uq 
polytrope, again excluding the small energies. 

Summarizing, we have shown that, for some classes of 
initial conditions, the QSS can be characterized by poly- 
tropic distributions. This appears to be the norm rather 
than the exception. This does not solve the problem of 
predicting the QSS structure as a function of the initial 
conditions since there is no theory to predict the poly- 
tropic index n, and there exist cases where the QSS is 
not even a polytrope (hence, polytropic distributions are 
not "universal attractors" ) . The problem of predicting the 
QSS would probably need a proper kinetic equation rather 
than a theory based on pseudo-entropy functionals (see 
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discussion in [2H])- However, we hope that these results 
can give some hint for the search of a more quantitative 
explanation. We have also shown that for other classes of 
initial conditions the QSS can be described by a Lynden- 
Bell distribution (see Fig. 36). The distinction between 



polytropic (Tsallis) and Lynden-Bell distributions may be 
related to the ergodic properties of the dynamics. When 
the system "mixes well" , the Lynden-Bell distribution is 
reached. Otherwise, the relaxation is incomplete and poly- 
tropic distributions may emerge. We conclude that both 
Lynden-Bell and Tsallis distributions may be relevant to 
describe QSSs depending on the efficiency of mixing as 
argued in ^TE\ . 
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